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FOREWORD 


TRACE  was  conceived  ».n  August  1961,  to  meet  the  Aerospace  Corporation 
requirements  in  the  fields  of  satellite  tracking  and  system  design.  Although 
under  continuing  development,  the  program  has  been  used  extensively  in  post¬ 
flight  orbit  determination  and  tracking  system  analysis. 

The  present  report  is  the  first  "complete"  description  of  the  program;  the 
previous  partial  descriptions  issued  in  March  and  December  1962  are  now 
obsolete  and  should  be  discarded.  The  information  herein  should  be  suffi¬ 
cient  for  most  users  of  the  program. 

TRACE  was  designed  by  M.  Bennett,  R.  J.  Mercer,  D.  Morrison, 

L.  Sachnoff,  and  C.  C.  Tonies;  the  principal  contributors  to  the  program 
include  the  originators  and  D.  A.  Adams,  C.  Christensen,  D.  Groves, 

K.  Hubbard,  S.  McDonald,  J.  Ostlie,  and  A.  Skulich. 
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ABSTRACT 


TRACE  is  a  multiple-purpose  satellite  orbit-determination  program  for  the 
IBM  7090  computer.  Its  applications  include:  (1)  prediction  -  the  generation 
of  a  satellite  trajectory  and  associated  ground  trace  and  station  sighting  data; 

(2)  orbit  determination  -  estimation  of  trajectory  parameters,  station  loca¬ 
tions,  and  observational  biases,  so  as  to  best  fit  a  set  of  observations;  and 

(3)  error  analysis  -  estimation  of  the  potential  accuracy  attainable  by  a 
tracking  system,  given  the  station  locations,  the  data  types,  rates  and 
quality,  the  uncertainties  in  the  model  parameters,  and  the  specifications 
of  the  nominal  orbit.  The  report  contains  the  objectives  of  the  program, 
some  theoretical  foundations,  the  equations  and  methods  employed,  the  struc¬ 
ture  of  the  program,  and  complete  instructions  for  its  use. 
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TRAJECTORY  GENERATION 


I.  Z 


Basic  to  all  applications  of  TRACE  is  the  space  vehicle  trajectory,  A  tra¬ 
jectory  is  defined  by  a  set  of  initial  conditions  together  with  the  differential 
equations  of  motion  which  reflect  the  earth's  gravitational  and  atmospheric 
forces  and  those  of  other  bodies  affecting  the  vehicle's  motion.  In  TRACE, 
the  trajectory  is  generated  in  an  inertial  rectangular  coordinate  system  by 
a  step-by-step  numerical  integration.  (Between  the  integration  points,  the 
trajectory  is  defined  by  an  interpolation  formula.  ) 

1.  Z.  1  Trajectory-Related  Output 

Trajectory-related  output  may  be  obtained  from  any  application  of  TRACE 
at  any  reasonable  set  of  time  points,  as  well  as  optionally  at  the  points  of 
equator  crossings,  apogee,  and  perigee.  Any  or  all  of  three  blocks  of  in¬ 
formation  may  be  selected  for  output  at  print  times.  They  are: 


1.  Z.  Z 


a.  Basic  trajectory  information  including  position  and 

velocity  components ,  spherical  coordinates,  ground 
trace,  and  altitude-.  * 

iSohe*-  r  '  ■  ■  *  VjW--, 

b.  Conic  section  elements'  computed  .from'  the  components 

of  position  and  velocity.^  V  *  * 

•  .  .  *  ••  ~ 

c.  Partial  derivatives  of  position  andsvelocity  components 
with' respect  to  initial  conditions:  and  dilferential  equation 
parameters.  (These,, quantiti’es';,a re  necessary  for  the 
other  applications,  of  TRACE  but  have  also^been  found 
useful  in  simple,  trajectory generation  where  effects  of 
initial  condition  or  parameter  errors  are  sought.) 

.  ' .  '■*'  xpiyl  '  y- i '  . 

Required  Input  ■'  ^  •  *’  • 


a.  Epoch  -  the  date  and  time  of  injection. 

b.  Initial  conditions  of  the  orbit.  Three  types  are  acceptable: 
(1)  inertial  rectangular  components  of  position  and  velocity, 
(Z)  spherical  coordinates  for  position,  together  with  the 
flight  path  angle,  azimuth,  and  magnitude  of  the  velocity 
vector,  and  (3)  elements  of  a  conic  section.  In  (Z),  either 
the  right  ascension  (inertial)  or  the  longitude  (referred  to 
Greenwich)  may  be  specified. 


c.  The  drag  parameter  CpA/W  and  a  choice  of  atmosphere 
model.  (The  ARDC  1959  model  is  used  unless  otherwise 
specified.  ) 

1.  2.  3  Optional  Input 

a.  Print  times  and  output  block  indicators.  {Trajectory- 
related  output  is  optional  for  any  of  the  applications  of 
TRACE;  in  the  simple  trajectory  generator  function,  this 
output  is  presumably  the  purpose  of  the  run.  ) 

b.  If  partial  derivatives  with  respect  to  certain  trajectory 
parameters  are  required,  these  parameters  must  be 
specified. 

c.  Many  model  constants  and  numerical  integration  param¬ 
eters  have  been  assigned  standard  values.  All  of  them 
may  be  changed  by  optional  input. 
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RADAR  DATA  GENERATION 


1.  3 

Simultaneous  with  trajectory  generation,  TRACE  can  produce  listings  of 
satellite  rise  and  set  times,  radar  coordinates,  and  many  related  quantities 
for  up  to  50  radar  stations,  provided,  of  course,  that  the  location  and  charac¬ 
teristics  of  the  stations  arc  supplied.  _  During  visibility  periods  (determined 
by  the  program),  the  desired  output  quantities 'are  computed  in  chronological 
order  and  stored  in  the  memory  until  capacity*  is-  reached,  at  which  point  the 
information  is  sorted  by  station  and  output.  The  process  is  repeated  as 
necessary  to  complete  the  listings.  Optionally  the  eight  quantities,  through 
Q  in  the  following  paragraph,  may  be  written  on  a  magnetic  tape  in  chron- 
ological  order  in  the  format  of  tracking  input  data. 

The  quantities  to  be  output  (in  the  station  listings  or  on  the  data  tape)  are 
selected  by  input,  and  include  range,  azimuth,  elevation,  range  rate, 
doppler  data  (P,  O,  P  and  Q) .  azimuth  rate,  elevation  rate,  range  accel¬ 
eration,  mutual  visibility  (up  to  eight  stations  only),  latitude,  longitude, 
surface  range,  altitude,  doppler  rate,  look  angle,  and  standard  deviations 
of  the  six  observational  quantities  R,  A,  E,  R,  A,  E. 

1.3.1  Required  Input 

a.  Station  location  data. 

b.  Control  information  for  each  station,  such  as  t.he  minimum 
and  maximum  elevation  angles,  maximum  range  of  visi¬ 
bility,  the  interval  (during  visibility  periods)  at  which 
computations  are  to  be  made,  and  the  start  and  stop  times 
for  visibility  testing  and  outpj.it. 

c.  The  list  of  output  quantities  desired  from  each  station. 


*  "Radar  station"  should  be  interpreted  here  generally  as  a  point  on  the 
surface  of  the  earth  associated  with  satellite  observations.  It  could  be 
a  camera  location,  or  the  location  of  a  point  observed  from  the  satellite. 

Optionally,  random  noise  may  be  added  to  the  same  eight  observational 
quantitie  s . 


1.  3.  2 


Optional  Input 


a.  Control  flags  to  indicate  that  only  rise  and  set  times  are 
to  be  generated,  and  that  the  generated  observational 
quantities  are  to  be  listed  chronologically  on  a  magnetic 
tape  in  the  format  of  input  data. 

b.  The  mean  and  standard  deviation  of  normally  distributed 
random  noise  to  be  added  to  the  generated  observational 
quantitie  s. 

c.  The  computed  value  of  elevation  is  altered  to  account  for 
atmospheric  refraction.  A  numerical  coefficient  may  be 
changed,  or  set  to  zero  if  no  correction  is  desired,  by 
optional  input. 

d.  Uncertainties  in  the  initial  conditions  of  a  trajectory  will 
be  reflected  in  uncertainties  in  generated  observational 
quantities.  A  vari ance -covariance  matrix  of  initial  con¬ 
ditions  must  be  input  if  the  computations  of  the  standard 
deviations  of  observational  quantities  are  selected. 
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1.4 


TRACKING 


Approximately,  TRACE  can  determine  the  trajectory  that  best  fits  a  set  of 
obse  rvations . 

1.4.  1  The  Tracking  Problem 

More  precisely,  the  trajectory  of  a  space  vehicle  depends  upon  the  initial 
conditions  of  the  motion  and  the  differential  equation  parameters  which 
appear  in  the  equations  of  motion.  From  the  trajectory,  one  may  compute 
at  the  observation  times  the  expected  values  of  the  recorded  observations. 
This  computation  further  depends  upon  the  locations  of  the  radar  stations 
and  biases  in  the  observations.  Thus,  the  computed  "observations"  are 
functions  of  parameters  of  four  types:  initial  condition,  differential  equation, 
station,  and  observation  parameters.  The  tracking  problem  is  to  solve  for 
the  set  of  parameters  that  minimizes  the  differences  between  the  computed 
and  measured  observations. 

Therefore,  TRACE  is  able  to  solve  for  such  quantities  as  the  ballistic  co¬ 
efficient  (a  differential  equation  parameter)  of  the  vehicle,  the  location  of  an 
observing  station,  and  the  presence  of  observational  or  time  biases  in  the 
.  data  reported  by  a  station,  in  addition  to  the  usual  initial  condition  param- 
eters.  {In  practice,  one  solves  only  for  a  selected  set  of  parameters  rather 

th'an'-all  possible  parameters.  ) 

.  ■ 

1 . 4 .  2  The  Tracking  Problem  Solution 

The  solution  is  an  iterative  process.  Initial  estimates  of  each  of  the  param¬ 
eters  must  be  provided.  Based  on  these  estimates,  the  "computed  obser¬ 
vations"  and  their  partial  derivatives  with  respect  to  the  parameters  are 
formed,  the  normal  matrix  is  accumulated,  and  measured  and  computed 
observations  are  differenced,  forming  the  "residuals."  The  residuals  are 
weighted  by  a  combined  scale  and  quality  factor,  checked  against  an  editing 
criterion,  and  the  sum  of  the  squares  of  the  weighted  residuals  is  accumu¬ 
lated.  When  all  of  the  observations  have  been  so  treated,  a  correction  to 


the  set  of  parameters  is  computed  and  applied,  and  the  process  is  repeated. 
The  root  mean  square  of  the  weighted  residuals  provides  the  measure  of 
convergence  of  the  process. 

The  solution  parameters,  which  are  derived  from  observations  containing 
random  errors,  must  be  regarded  as  estimates  of  the  true  parameters. 
Under  certain  conditions  (that  the  observational  model  is  correct,  that  the 

observational  errors  are  independently  distributed  with  mean  zero  and 
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variance  cr  ,  and  that  the  weighting  factor  used  is  cr  )  the  inverse  of  the 

normal  matrix  is  the  variance-covariance  matrix  of  the  parameters.  Thus 

the  solution  process  provides  an  estimate  of  the  uncertainties  in  the  derived 

parameters. 

Two  types  of  conditions  may  be  imposed  upon  the  solution  of  the  minimiza¬ 
tion  problem:  bounds  upon  the  magnitude  of  the  computed  corrections  may 
be  given,  and  linear  constraints  among  the  corrections  may  be  specified. 

The  former  is  used  to  assure  covergence,  and  the  latter  may  represent 
physical  requirements,  such  as  the  fact  that  the  difference  in  the  latitudes 
of  two  stations  is  accurately  known. 

The  output  includes  rhe  rms  of  the  residuals  (optionally  the  residuals  may 
also  be  reported  by  station)  for  the  current  iteration,  the  current  and  cor¬ 
rected  values  of  the  parameters,  the  rms  residual  that  is  predicted  for  the 
next  iteration,  and  the  standard  deviations  of,  and  the  correlations  among, 
the  parameters  (obtained  from  the  variance -covariance  matrix). 

Further  optional  output  includes  the  individual  residuals,  partial  derivatives 
of  observations  with  respect  to  parameters,  and  trajectory  information  at 
observation  times. 


1.4.  3 


Required  Input 


a.  The  list  of  initial  condition,  differential  equation,  station 
location,  and  observational  parameters  for  which  the  pro¬ 
gram  is  to  solve,  an  initial  estimate  of  each,  and  a  bound 
(if  necessary)  upon  the  magnitude  of  correction  that  is  to 
be  permitted 

b.  Locations  of  all  observation  stations 

c.  The  observational  data.  Many  types  of  data  are  acceptable, 
but  eight  of  these  (range,  azimuth,  elevation,  range  rate, 
and  the  four  doppler  quantities,  P,  Q,  P,  Q)  are  regarded 
as  basic;  the  other  data  types  are  first  converted  to  the 
above  set. 

d.  Weighting  factors  for  the  basic  data  types  for  each  observ¬ 
ing  station 

1.4.4  Optional  Input 

a.  The  maximum  number  of  iterations  in  the  differential 
correction  process  may  be  specified. 

b.  A  refraction  correction  is  applied  to  elevation  observations. 
A  coefficient  in  this  correction  may  be  modified,  or  set  to 
zero,  by  optional  input. 

c.  The  names  of  up  to  nine  stations  for  which  residuals  are  to 
be  reported  may  be  input. 

d.  If  linear  constraints  among  the  parameters  exist,  a  con¬ 
straint  matrix  must  be  input. 

e.  The  level  of  residuals  above  which  data  points  are  to  be 
discarded  may  be  input. 
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STATISTICAL  ANALYSIS 


1.  5 

As  already  mentioned,  under  certain  assumptions  the  inverse  of  the  normal 
matrix  contains  information  as  to  the  uncertainties  with  which  parameters 
are  determined  by  a  tracking  system.  This  aspect  of  orbit  determination  is 
exploited  in  TRACE  to  provide  its  system  analysis  capability.  (Note  that  the 
normal  matrix  involves  only  partial  derivatives,  not  residuals,  and  thus  the 
performance  of  a  system  can  be  analyzed  without  recourse  to  actual 
observations.  ) 

1.5.1  Assumptions 

The  assumptions  are  that  observational  errors  are  independently  distributed 

2  - 1 

with  mean  zero  and  variance  <r  ,  and  that  cr  is  used  as  the  weighting  factor 
in  forming  the  normal  matrix.  Under  these  conditions,  the  inverse  of  the 
normal  matrix  is  a  variance-covariance  matrix  of  the  parameters  being 
estimated, 

Insofar  as  these  parameters  are  differential  equation  or  radar  parameters, 
their  variances  and  covariances  satisfactorily  describe  their  uncertainties. 
However,  the  uncertainties  in  the  motion  of  a  vehicle  are  not  adequately 
described  by  the  variances  and  covariances  of  the  initial  conditions  and 
differential  equation  parameters  of  the  trajectory.  Rather,  trajectory  un¬ 
certainties  are  better  reported  in  terms  of  orbit  plane  coordinates,  or  conic 
section  elements,  or  such  related  quantities  as  period,  apogee,  and  perigee 
distance.  Cartesian  and  spherical  coordinate  variance -covariance  matrices 
are  also  available. 

A  further  sophistication  arises  from  the  assumption  that  the  values  of  some 
of  the  parameters  used  in  the  calculations,  but  not  being  estimated  by  the 
differential  correction  process,  are  also  uncertain,  thereby  inducing  un¬ 
certainties  in  the  differentially  corrected  parameters  and  in  the  trajectory. 
(This  is  a  very  common  situation;  most  tracking  programs  do  not  solve  for 
basic  constants  and  station  locations,  but  their  current  values  must  be  some¬ 
what  uncertain.)  Such  parameters  are  referred  to  as  "Q-parameters"  in 
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distiction  to  "P-parameters,  "  which  are  those  being  estimated  by  differential 
correction.  TRACE  will  simultaneously  report  P-parameter  and  trajectory 
uncertainties  with  and  without  Q-parameter  effects.  The  matrix  of  derivatives 
of  the  P-parameters  with  respect  to  the  Q-parameters  can  also  be  output. 

1.5.2  Mechanics  of  Application 

The  mechanics  of  this  application  of  TRACE  are  as  follows:  as  the  trajectory 
is  generated,  the  program  determines  (and  outputs)  periods  of  visibility'  from 
each  station.  At  the  prescribed  interval,  while  the  vehicle  is  visible,  the 
partial  derivatives  of  the  observations  with  respect  to  the  P-  and  Q-parameters 
are  computed,  weighted,  and  accumulated  (in  double  precision)  into  the  normal 
matrix.  At  the  specified  output  points,  the  desired  covariance  matrices  are 
computed  and  output,  along  with  such  trajectory- related  quantities  as  may 
have  been  selected. 

1.5.3  Required  Input 


a.  The  list  of  up  to  30  parameters,  of  which  up  to  15  may  be 
trajectory  parameters 

b.  The  set  of  output  times  and  the  list  of  output  variance - 
covariance  matrices  desired.  Optionally,  only  the  standard 
deviations  (square  root  of  the  diagonal  elements)  can  be 
printed. 

c.  Station  location  information 

d.  Control  information,  for  each  station,  such  as  the  minimum 
and  maximum  elevation  angles,  maximum  range  of  visibility, 
the  interval  (during  visibility  periods)  at  which  computa¬ 
tions  are  to  be  made,  and  the  start  and  stop  times  for 
visibility  testing  and  output 

e.  The  list  of  data  types  reported  by  each  station.  Eleven 
types  are  possible:  range,  azimuth,  elevation,  range  rate, 
the  four  doppler  quantities  (P,  Q,  P,  Q) ,  argument  of  the 
latitude,  the  orthogonal  angle  measured  from  the  equator 
to  the  position  of  the  vehicle  in  the  plane  containing  the 
radius  vector  and  the  vector  normal  to  the  orbit  plane,  and 
geocentric  distance.  (The  last  may  be  used  to  simulate 
altitude  observations,  as  height  and  geocentric  distance 
have  a  nearly  constant  difference.  ) 

f.  Standard  deviations  for  each  type  of  data 
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a.  If  Q-parameter  effects  are  desired,  each  parameter  must 
be  designated  as  being  of  type  P  or  Q. 

b.  An  input  covariance  matrix  is  required  for  the  set  of 
Q-paramete  rs . 
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SECTION  2 


THEORY 


2.  1  INTRODUCTION 

Section  1  contains  general  descriptions  of  the  various  applications  of  TRACE 
to  orbit  determination  problems.  In  Section  2,  more  precise  mathematical 
statements  of  these  problems  are  provided  and  the  capabilities  of  TRACE 
are  discussed.  The  emphasis  in  this  section,  however,  is  upon  theoretical 
aspects  and  functional  relations;  the  particular  equations  and  methods  used  in 
the  program  are  set  forth  in  Sect:'''  7  Again,  the  applications  are  treated 
in  the  order  of  increasing  scope,  but  not  with  uniform  thoroughness.  Topics 
that  are  possibly  less  familiar  have  been  emphasized,  whereas  more  familiar 
problems,  such  as  the  numerical  solution  of  differential  equations,  have  been 
largely  ignored. 
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THE  TRAJECTORY  AND  ITS  PARTIAL  DERIVATIVES 


2.  2 

The  trajectory  of  a  space  vehicle  is  defined  by  the  (differential)  equation 
of  motion 

X  =  -^  +  F 


together  with  the  initial  values  X(t  )  =  X  and  X(t  )  =  X  .  Here.  X  is  a 

o  o  o  o 

3 -vector  of  rectangular  components  (x,  y,  and  z)  of  position  in  an  inertial 

coordinate  system,  a  dot  represents  a  time  derivative,  r  =  X  - 
2  2  2  1  /  2 

(x  +  y  +  z  )  ,  p.  is  the  gravitational  constant  (GM)  of  the  earth,  and  F 

(a  vector)  represents  the  perturbing  accelerations  upon  the  vehicle. 


One  application  of  TRACE  is  merely  to  solve  this  differential  equation. 

solution  X(t),  X(t)  is  generated  numerically  at  time  points  t  =  t.  (j  =  0, 

J 

2,  .  .  .  )  and  defined'  at  t  ^  t.  by  an  interpolation  formula. 


The 

1, 


The  more  sophisticated  applications  of  TRACE  require  the  sensitivity  (as 
expi  3ssed  by  partial  derivatives)  of  the  trajectory  to  its  initial  conditions 
and  other  parameters. 


Obviously  X  is  a  function  of  p.  which  is  an  example  of  a  ’differential  equation 
parameter."  Other  such  parameters  (ballistic  coefficients,  oblateness 
coefficients,  etc.  )  may  appear  in  F.  Furthermore,  the  solution  depends  on 
the  initial  conditions  X^  and  Xq,  which  in  turn  may  be  computed  from  "initial 
condition  parameters.  "  If  we  let  vectors  of  these  typos  of  parameters  be 
represented  by  (3  and  a  respectively,  we  may  show  the  functional  relations 
in  Eq.  ( 1 )  as 


X  =  -  Hy  +  F(X,  X,  P,  t) 
r 


(2) 
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or 


X  =  X{X,X,p,t) 


with 


X(tQ)  =  Xq(q)  ,  X(tQ)  =  Xo(q)  (2a) 

(F  and  X  will  be  functions  of  X  whenever  drag  forces  are  present.)  The 
dependence  of  the  solution  upon  the  parameters  can  be  indicated  as 

X(t)  =  X(a ,  (3 ,  tQ,  t) 

and  similarly  for  X(t).  (In  fact,  the  solution  can  be  given  by  the  integral 
equations 


X(a,p,to,t)  =  Xq(q)  +  j^X[X( a,p,to,t"),  X(a,p,to,t"),P,t"]  dt" 


and 


X(a,p,to,t)  =  XQ(a)  +  C  X(a ,  (3 ,  t  ,  t'  )  dt' 

J  o 


=  X  (a)  +  (t  -  t  )X  (a) 

O  o  o 

•t  rt 


f  f  X[X(a,  (3,  tQ,  t"),  X(a,p,to,t"),  P  ,  t" ]  dt"dt' 
J  o 


X(t)  =  XQ(a)  +  (t  -  tQ)Xo(a) 


+  J^  (t  -  t"  )X[X(a,pfto,  t"),  X(a,p,to>  t"),  P,t"]dt' 


(3) 
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which  are  hardly  suitable  for  computations,  but  which  do  show  the  functional 
relations  more  explicitly.  ) 


3X  ax  ax 

We  are  now  in  a  position  to  show  how  partial  derivatives  7^-,  and  7^—  , 

which  measure  the  sensitivity  (to  first  order)  of  solutions  to  variations  ° 
in  the  trajectory  parameters  a,  (3,  and  tQ,  are  obtained.  (These  partial 
derivatives  are  extensively  used  in  other  applications  of  TRACE,  but  are 
also  often  of  interest  in  their  own  right.  ) 

If  we  differentiate  Eqs.  (2)  and  (2a)  with  respect  to  a,  interchange  orders 

3X 

of  differentiation,  and  use  the  notation  X  for  rr — ,  we  obtain 

a  9a 


ax  px\  9f]9X  ,  aF  ax 

3a  ^3Xy  ^3  J  3X  3a  3X  3a 


or 


X  -  r  3  (  M-X\  8F 
Xa  -  ,L3X^7-J+  3X 


X  +  £?  X 
a  3X  a 


(4) 


3X  3X 

with  initial  conditions  X  (t  )  =  — ^ —  and  X  (t  )  =—7; — .  (See  Paragraph  2.2.1). 

a  o  3a  a  o  3a  r  ' 

Equation  (4)  is  called  a  "variational  equation.  "  It  is  obviously  a  second-order 
linear  vector  differential  equation  whose  solution  is  the  vector  of  partial 

3X 

derivatives  X  =  — —  of  the  components  of  position  with  respect  to  the  initial 
a  3a  r  •  3X 

condition  parameter  a.  In  the  course  of  solving  Eq.  (4),  X^  =  -7^-  will  also 

be  obtained.  Such  an  equation  can  be  derived  for  each  initial  condition 

parameter. 

One  can  also  obtain  Eq.  (4)  by  differentiating  the  integral  Eq.  (3)  with  respect 
to  a , 


X  =  X  (t  )  +  (t 
a  av  o 


-  t  )X  (t  )  +  f  (t  -  t"  )(|$  +  ^  |^)  dt" 

o  a  o  J  t^  \3X  3a  gX  3a/ 


(5) 
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and  noting  that  Eq.  (5)  corresponds  to  Eq.  (4)  in  exactly  the  same  way  that 
Eq.  (3)  corresponds  to  Eq.  (2).  Reference  1,  which  initiated  the  authors' 
use  of  variational  equations,  follows  the  integral  formulation,  but  is  more 
concerned  with  interplanetary  applications. 

The  variational  equations  for  initial  time  tQ  are  of  the  same  form,  but  with 
different  initial  conditions 


3  F 

ax 


X. 


-I  o 


+ 


^x 

ax  t 


{t 


-Xo(a) 


Xt  <‘o> 
o 


-X  {a' 
o 


These  are  derived  by  differentiating  the  integral  equation  (Eq.  3),  which  best 
shows  the  dependence  upon  tQ,  with  respect  to  tQ. 

The  variational  equations  for  a  differential  equation  parameter  (3  are 


3  ax 

r  > 


X  +  X  + 

P  ax  P  as 


W  =  VV  =  0 


(6) 


As  a  source  of  partial  derivatives,  variational  equations  give  results  that  are 

more  accurate  than  analytic  derivatives  (which  assume  two-body  motion),  and 

are  more  rapidly  generated  than  difference  quotient  approximations.  The 

1"  3  /  pX  - __1 

greater  speed  derives  from  the  fact  that  the  terms 


a 

ax 


+ 


a  f 


ax 


and 


of  Eq.  (4)  are  identical  in  all  the  variational  equations;  only  the  non- 

3  F 

homogeneous  term  —  of  Eq.  (6)  varies  with  the  particular  parameter. 


3F 

ax 


A  further  advantage  of  the  variational  equations  is  that  they  permit  the  use  of 
the  difference  quotient  technique  as  a  checking  device.  The  two  methods 
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must  produce  partial  derivative  estimates  that  are  in  substantia]  agreement; 
the  lack  thereof  would  indicate  the  presence  of  a  blunder.  While  the  test  is 
hardly  foolproof,  it  is  valuable  and  should  rot  be  overlooked. 


2.  2.  1  Derivative  with  Respect  to  a  Vector 

The  indication  of  a  derivative  with  respect  to  a  vector  is  a  very  convenient 
notational  device  for  representing  partial  derivative  matrices  and  chain 
rule  differentiation.  The  following  conventions  are  observed  throughout: 

a.  A  “vector"  is  a  column  vector;  a  row  vector  will  be 
described  as  such  or  denoted  as  a  transposed  vector. 
Example:  (x,  y,  z)  =  X^1. 

b.  The  derivative  of  a  vector  with  respect  to  a  scalar  is  a 
vector. 

c.  The  derivative  of  a  scalar  with  respect  to  a  vector  is  a 
row  vector. 

d.  The  derivative  of  a  vector  with  respect  to  a  vector  is  a 
matrix. 

Example:  If  F  is  a  vector  function  of  a  vector  variable  X, 

rj  jp  £ 

then  -r-v  is  the  matrix  of  partial  derivatives  whose  i  -  j 
3X  OFi 
element  is  . 

a  Aj 

Note  the  neatness  of  the  following  example:  Suppose  X(t)  is  the  vector 

T 

(Xj(t),  x^t),  .  .  .  ,  xn(t)]  and  y  is  a  scalar  function  of  X;  y  =  f[xJ(t), 
x2(t) . xn(t)]  =  f[X(t)].  Then 


-  JLL  +  .11  4.  at  dxn  _  9f  dx 

dt  3x.  dt  9x_  dt  '  3x  dt  ~3X  dt 

1  2  n 


Here  is  by  convention  a  row  vector,  a  column  vector,  and  the 

juxtaposition  of  the  two  indicates  the  desired  scalar  product. 
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BASIC  ORBIT  DETERMINATION 


2.  3 

The  basic  orbit  determination  problem,  as  outlined  in  Section  1.4,  is  that  of 
finding  values  for  a  set  of  parameters  from  an  observational  model  so  as  to 
minimize,  in  the  sense  of  weighted  least  squares,  the  differences  between  the 
measured  observations  and  the  corresponding  quantities  computed  from  the 
model. 

The  model,  as  constituted  in  TRACE,  includes  the  trajectory  of  the  vehicle 
(and  thus  the  initial  condition  and  differential  equation  parameters),  the 
locations  of  the  observing  stations,  and  constant  bias  errors  in  their  instru¬ 
ments  or  their  clocks.  In  practice,  one  determines  values  for  only  a  selected 
subset,  P,  of  the  parameters  of  the  model. 

The  weighting  factors  are  necessary  to  assign  the  proper  relative  significance 
to  observations  of  different  types  and  quality. 

The  basic  orbit  determination  problem  may  now  be  restated:  Given  a  set  of 
n  normalized  observations  (multiplied  by  an  appropriate  weighting  factor', 
which  are  collectively  denoted  by  the  n-vector  (m  for  "measured''),  and 

a  model  from  which  the  corresponding  (similarly  weighted)  quantities  O  can 
be  computed  as  functions  of  parameters  P,  determine  values  of  P  so  that 
II  O  -  O  (  P)  1 1  ^  is  minimized. 

Suppose  that  an  approximate  solution  P^  is  known.  (Approximate  initial 

conditions  will  be  available  either  from  design  information  or  preliminary 

orbit  determination  methods.  )  We  expand  O  (P)  in  a  Taylor  scries  to  first 

order  about  P  and  obtain 
o 

l|Om  -  Oc(P)||2  =  |IOm  -  Oc(Po)  -  A-  AP||2  (8) 

DO 

c  . 

to  be  minimized,  where  A  =  ■  — p-  is  a  matrix  of  normalized  partial  derivatives 
evaluated  at  P  =  Pq.  The  partial  derivatives,  with  respect  to  trajectory 


a°c  90  x 

parameters,  are  computed  from  the  chain  rule  formula  ^ 

^  ^  C/  9  **  9  x 

where  is  the  matrix  of  solutions  to  the  variational  equations.  The  matrix 

aoc  9P  aoc 

and  those  columns  of  that  represent  derivatives  with  respect  to 

station  parameters  are  computed  directly  from  geometrical  relations. 


The  differences,  O  (P  )  =  O  -  O  (P  ),  between  the  normalized  observa- 

mc  o  m  c  o 

tions  and  the  corresponding  quantities  computed  from  assumed  values  Pq 
are  called  "residuals";  they  will  be  due  to  the  presence  of  random  observa¬ 
tional  errors,  inadequacies  in  the  model,  and  incorrect  values  for  the  model 
parameters. 


The  above  statement  of  the  weighted  least  squares  (WLS)  problem  conceals 

the  weighting  factors,  which  could  have  been  explicitly  included  in  the 

formulation  as  the  elements  of  a  diagonal  matrix  {This  notation  is 

l/2  T  l/2 

chosen  in  order  to  simplify  subsequent  equations  in  which  (W  )  W  =  W 
appears  frequently.  )  Then  the  quantity  to  be  minimized  would  have  been 
I  W  (°mc  -  A-AP)  I  ,  with  Omc  and  A  representing  actual,  not  weighted, 
values.  Since  TRACE  is  restricted  to  independent  observations  for  which 
W  is  diagonal,  we  have  chosen  {for  the  sake  of  simplicity,  not  deception)  to 
include  the  weighting  factors  with  the  elements  of  O  and  A.  In  Section 
2.  5  on  statistical  aspects,  the  weighting  matrix  is  discussed  and  displayed. 


It  should  be  noted  that  the  solution  of  the  WLS  problem  docs  not  produce  "true''1 
values  for  model  parameters  —  it  produces  only  best  fit1’  values.  Any 
further  conclusions  of  a  statistical  nature  regarding  the  WLS  solution  require 
additional  assumptions  regarding  the  model  and  the  character  of  the  random 
observational  errors.  These  topics  will  be  discussed  later. 


Inasmuch  as  the  original  nonlinear  WLS  problem  has  been  replaced  by  an 

approximate  linear  problem  (that  of  finding  a  correction  sector  AP  so  that 

O  { P  )  -  A-AP  ^  is  minimized),  we  must  not  expect  that  P  ~  P  +  AP 
m  c  o  o 

will  be  a  solution  of  the  original  problem;  rather,  an  iterative  process  is 

indicated.  O  (P  )  measures  the  degree  to  which  an  orbit,  computed 
me  o  1 
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which  must  be  zero  if  AP  minimizes  f(AP). 

Proof  2  - 
Let 

A^A  *  AP  =  A^O  .  Then  for  any  AP'  *  AP 
me 


f{AP'  )  =  t  I  A  -  AP  -  O  +  A(AP'  -  AP) 

'  '  11  me 

=  ||  A*  AP  -  O  ||Z  +  2[A  •  AP  -  °mc]TtA^Ap/  "  AP)^ 
ii  me 11  mc 

+  ||A(AP'  -  AP)||2 

=  £(  AP)  +  Z(ATA  •  AP  -  ATOmc)T{  AP'  -  AP) 

+  ||A(AP'  -  AP)  |  |  2 

=  £(AP)  +  |  |  A(  AP'  -  AP)  |  |  2 
>  f(  AP) 

from  which  we  see  that  AP- minimizes  f(AP). 


2-10 


CONSTRAINED  AND  BOUNDED  LEAST 


2.  4 

SQUARES  SOLUTIONS 

Two  distinct  types  of  restrictions  upon  the  solution  of  the  weighted  least 
squares  problem  may  be  necessary  or  desirable.  Constraints  among  the 
parameters  may  be  a  part  of  the  physical  problem,  and  bounds  upon  the 
magnitude  of  the  corrections  may  be  computationally  desirable. 

2.  4.  1  Constraints 

An  example  of  a  physical  constraint  among  parameters  would  be  precise 
knowledge  of  the  relative  locations  of  two  nearby  radar  stations.  If  their 
actual  locations  were  among  the  parameters  P  in  a  differential  correction, 
it  would  be  important  to  constrain  the  corrections  AP  so  that  the  radar 
relative  locations  were  preserved.  This  is  accomplished  in  TRACE  by  intro¬ 
ducing  linear  constraints  in  the  form 


AP  =  B  •  AP'  f  C 


(10) 


where  B  is  a  rectangular  matrix  and  AP'  is  a  reduced  set  of  independent 
parameters;  by  solving  the  WLS  problem  in  terms  of  AP'  ;  and  by  using 
Eq.  (10)  to  obtain  the  constrained  corrections  AP.  Solving  the  WLS  problem 
in  terms  of  AP'  requires  the  minimization  of  j  I  A  •  AP  -  O  11^  subject 
to  the  constraint  Eq.  (10),  or  therefore  the  minimization  of  |jA-(B’AP'  +  C) 
-  Orri(  i  -  i|(AB)'  AP'  -  (O^  +  AC)j  j  .  The  solution  of  the  linear  system 

(AB)T(AB)AP'  =  ( AB)T(  +  AC)  (11) 


gives  the  required  minimum. 

2.4.2  Bounds 

Under  fairly  common  conditions,  such  as  inadequacies  in  the  observational 
model  or  a  poor  initial  approximation  Pq,  the  observed  !l^mcll  wiH  fail  to 

Radar  station  may  refer  to  any  point  on  the  surface  of  the  earth  associated 
with  an  observation. 
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follow  the  predicted  jj  CF^jj  or  may  even  diverge.  In  the  presence  of  such 
manifestations  of  nonlinearity,  it  may  be  necessary,  in  order  to  assure 
eventual  convergence,  to  solve  the  WLS  problem  at  each  iteration  with  a  side 
condition  bounding  the  magnitude  of  the  correction  vector  AP. 

If  we  refer  to  the  reciprocals  of  the  bounds  g.  collectively  as  the  diagonal 
matrix  G,  the  restricted  problem  becomes  that  of  minimizing  ]  A  •  AP  - 
subject  to  the  bounding  condition  |  j  G  •  AP  <  1.  (If  G  •  AP  |  j  ^ 

<  1,  then  for  each  component  j  Ap^  |  <  g..  ) 

The  constraint  has,  in  a  two-parameter  example,  a  simple  geometrical 
description.  The  constrained  problem  is  to  find  a  minimum,  over  all  AP 
within  the  ellipse  defined  by  g  ^  and  g^.  of  the  surface  i(AP)  =  [  J  A  •  AP  - 
(See  Figure  1 .  ) 


Figure  1.  Two-Parameter  Constraint  Ellipse 


{An  elliptic  rather  than  circular  region  is  used  to  account  for  the  range  of 
magnitudes  of  the  various  parameters.  ) 

If  the  unconstrained  solution  is  not  within  the  ellipse  ||G-AP]|^  =  1,  then  we 
invent  a  new  function  F(  AP)  =  f(  AP)  +  z  |  G  •  AP  \  J  ^  to  be  minimized.  The 
minimum  point  AP'(z)  is  found  as  the  solution  of 

(ATA  +  z.GTG)AP  -  AT0  (12) 

me 

As  z  increases,  the  minimization  of  F  will  require  smaller  and  smaller 
values  AP'  (z).  (More  precisely,  it  will  be  shown  that  ||G*AP'(z)||  is  a 
decreasing  function  of  z.  )  in  particular  we  can  find,  by  a  search  and  inter¬ 
polation  procedure,  a  value  z'  of  z  such  that  I)G*AP'(z/)||^=1*  But  for  z'  , 
the  minimization  of  F  =  f  +  z'  |  j  C-  •  AP'  (z '  )  1 1  =  f  +  z'  is  equivalent  to  mini¬ 
mizing  f,  since  they  differ  only  by  the  constant  z'  .  Thus  we  have  found  the 
point  AP'(z'  )  which  minimizes  f(AP)  along  the  bounding  ellipse. 

We  will  also  show  that  f(AP'(z))  is  an  increasing  function  of  z.  Thus,  any 
interior  point  of  the  ellipse  corresponds  to  larger  values  of  z  and  of  f,  and 
therefore  the  constrained  minimum  point  is  on  the  boundary  and  is  the 
solution  AP'(z')  of  (A  J  A  +  z'GTG)AP  =  ATOrnc  for  which  ]  |  G  *  AP'  ( z'  )  1 1  2  =  1 . 

The  monotonic  decreasing  character  of  '  |  G  •  AP  |  j  as  a  function  of  z  is  shown 
as  follows.  (Primes  have  been  dropped  throughout  these  proofs.  )  If  we 
diffe  rentiate 


(ATA  +  zGTG)AP(z)  “  Al°mr 


(13) 


with  respect  to  z  we  obtain 


(ATA  +  zGTG)^(AP)  +  (GTG)AP  =  0 


(14) 
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-f-(AP)  =  -(A  A  +  zG  G)  *(G  G)Al 


This  expression  is  needed  in  the  equation  for  ^-|  |  G  •  AP  1 1  ^ ,  which  will  be 
shown  to  be  negative  for  all  z  >  0.  We  have 


A||G-  AP|i2  =  2(APT|(GTG)(i  AP) 


=  -2APT(G1G)(A1A  +  zGTG)‘ l(GfG)AP 


T  T  - 1 

Since  (A  A  +  zG  G)  is  positive  definite  for  positive  z> 


A||G' ap||2<  o 


whenever  (G  G)  •  AP  *  0. 

The  monotonic  increasing  character  of  f[AP(z)]  is  similarly  established  by 
showing  that.  >0,  as  follows: 


3f  d  ( A  P ) 

AP]  dz 


=  [z|ATAAP  -  ATOiric)J,]  ^-  ( A  ^  A  +  zGTG)"1(G1G)AP 

=  -2|"(ATA  +  zGTG)AP  -  AT0  -  ztG^OAp]1  [(A1  A  +  zGTG)'‘iGTG)AF 


But  since  (ATA  +  zGTG)AP  =  AfO 

me 


^  =  +2zAPT(GTG)(ATA  +  zGTG)“1(GTG)AP 


(19) 


J  /*  rT% 

and  >  0,  whenever  (G  G)AP  ±  0. 

2.  4.  3  Solution  of  the  Linear  System 

T  T  T 

The  solution  of  the  linear  system  (A  A  +  zG  G)AP  =  A  O  (and  the  inversion 

rj.  rp  ID  C 

of  the  coefficient  matrix  C  =  A  A  +  zG  G)  is  accomplished  by  a  special 
method  akin  to  that  known  classically  as  the  square  root  method.  (See 
Reference  2.  )  It  is  a  finite  (noniterative)  method,  applicable  only  to  sym¬ 
metric  matrices,  and  is  based  on  the  fact  that  a  symmetric  matrix  can  be 

T 

decomposed  as  a  product  of  the  form  C  =  LDL  where  L  is  a  lower  triangular 
matrix  with  (-  1)  as  diagonal  elements,  and  D  is  a  diagonal  matrix. 

In  such  a  representation  det  (L)  =  ±1  and  det  ( D)  =  det  (C).  Therefore  L  * 
exists  (and  also  has  the  form  of  L)  and  D  has  no  zero  elements  if  C  is  non¬ 
singular.  Therefore  two  equivalent  forms  are 

(1)  L"lC(LT)‘i  =  L"1C(L‘1)T  =  D 

(2)  C'1  =  (L“1)TD“1L"1 
or 

( 1 7  )  SCST  =  D 

(27  )  C“l  =  STD'lS  where  S  =  L~  1 

Thus  we  see  that  the  inversion  of  C  and  the  solution  AP'  =  (C  )A  O 

-j-  rn  c 

require  matrices  S  and  D  such  that  (!'  )  SCS  =  D. 
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A  bordering  technique  is  used  to  find  S  and  D.  At  the  stage  suppose 
that  the  k^1  order  principal  minors  of  S  and  D  have  been  found.  The  (k  +  l)St 
order  minors  require  the  vector  W  and  the  scalar  b  so  that 


d\  /S, 


-l/\d 


W\  /D 


c/  \0 


th  /Ck  d\  s 

where  C-  is  the  k  and  I  ^  j  the  (k  +  1)  order  minors  of  C.  It  is 

\d  a/ 

easily  verified  that  the  required  W  and  b  are 


w  =  STD->skd 


b  =  a  -  W  d 
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THE  STATISTICS  OF  ORBIT  DETERMINATION 


2.  5 

In  the  process  of  orbit  determination  by  the  method  of  weighted  least  squares 
(WLS),  no  assumptions  regarding  the  statistics  of  the  observational  errors 
need  be  made.  In  this  case,  no  statistical  conclusions  can  be  drawn  from 
the  results,  and  the  justifications  of  the  method  are  simply  that  it  minimizes 
residuals  (in  the  sense  of  WLS),  and  that  it  works  in  practice. 

On  the  other  hand,  if  two  common  assumptions  are  made,  namely  (a)  that  the 

observational  errors  are  independent  with  mean  zero  and  known  variance 

<r^,  and  (b)  that  the  multiplicative  weighting  factor  associated  with  each 
1  _  \ 

observation  is  a  ^  ,  then  the  inverse  normal  matrix  is  a  variance -c ovariance 
matrix  (often  abbreviated  "covariance  matrix")  of  the  parameters  being 
determined.  This  matrix  depends  only  on  the  partial  derivatives  of  the 
observations  with  respect  to  the  parameters,  thus  permitting  statistical  analysis 
of  a  tracking  network  in  the  absence  of  actual  or  simulated  observations. 

The  details  are  covered  in  the  next  few  sections.  The  relation  of  WLS  orbit 
determination,  as  in  TRACE,  to  other  criteria  (minimum  variance  and 
maximum  likelihood)  is  also  discussed. 

2.  5.  1  The  Variance-Covariance  Matrix 

We  assume  that  the  vector  of  measured  observations  O  is  the  true  value 

m 

O  (PJ  plus  a  random  error  «.  Our  linear  approximation  to  O^(P)  is 

O  (P)  =  O  (P  )  +  A-  AP  (22) 

C  CO 


and  the  residual  vector  is 


O  =  O  -  O  (P  ) 
me  m  co 


A  •  APt  +  e 


(23) 
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whe  re 


First,  we  show  that  AP'  is  an  unbiased  estimate  of  the  true  value  AP^.  (By 
this  we  mean  that  although  AP'  is  a  random  quantity  since  it  depends  upon 
the  residuals,  and  thus  upon  the  observational  errors,  the  expected  value  of 


AP'  is  the  true  value  AP^.  ) 


AP'  =  (ATWA)'1ATWO 


=  (ATWA)  1  ATW(A  •  APt  +  0 


=  APfc  +  (ATWA)'lATW( 


E(AP')  =  APt  +  E[(ATWA)'1ATW(]  =  APt 


by  the  linearity  of  E(  *  )  and  on  the  assumption  that  E(«)  =  0. 

The  vector  6P'  =  AP'  -  AP^  would  be  the  deviation,  due  to  random  errors, 
of  the  solution  AP'  from  the  true  value  AP^;  it  has  been  shown  to  have  the 
expected  value  zero. 

V/hat  is  now  the  expected  value  of  the  square  of  the  deviations  or  of  the  product 

T 

of  two  components  thereof?  The  answers  are  summarized  in  E(6P'  6P'  ), 

which  is  by  definition  the  estimated  covariance  matrix  C(P' )  of  the 


parameter  s* 


E(6P'  6P'  T)  =  E[(ATWA)  1ATW(mT)WA(ATWA)'1]  ,  (27) 


in  which  we  used  the  symmetry  of  the  matrices  W  and  A  WA,  or 


C(P')  =  (A  1  WA)  1 ATW  Z  WA(ATWA)‘  1 


This  is  the  general  form  of  the  covariance  matrix  for  a  WES  estimate  of  the 
parameters.  If,  however,  Z  is  diagonal,  as  per  the  first  assumption,  making 
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(29) 


it  possible  to  choose  W  =  Z  \  as  per  the  second  assumption,  to  be  the 
diagonal  weighting  matrix,  then  the  very  great  simplification 

C(P')  =  (ATWA)~l  =  (A1  Z^Af  1 

is  the  result.  This  is  the  basic  covariance  matrix  calculated  in  TRACE. 

2.  5.  2  Minimum  Variance  and  Maximum  Likelihood  Estimates 

In  most  instances  of  orbit  determination  from  observations,  the  method  of 
weighted  least  squares  (WLS)  requires  no  statistical  justification;  indeed 
there  is  none.  Its  aim  is  simply  to  produce  fits  and  predictions  of  acceptable 
quality.  In  other  applications,  such  as  systems  analysis  and  design,  statis¬ 
tical  conclusions  are  sought.  These  are  commonly  based  on  minimum 
variance  (MV)  or  maximum  likelihood  (ML)  estimations.  (MV  is  also  called 
"Markov.11)  The  purpose  of  this  section  is  to  describe  in  general  terms  the 
assumptions  governing  MV  and  ML  techniques  and  to  relate  them  to  the 
basically  WLS  method  in  TRACE. 

The  MV,  or  Markov,  estimate  of  AP  is  that  linear  unbiased  estimate  that 
minimizes  the  diagonal  terms  {the  variances)  of  the  variance -covariance 
matrix  of  parameters  P.  (See  Reference  3.)  The  formulas  are 

AP.  ,  .r  =  { AT  Z~1A)~^AT  Z_10  (30) 

MV  '  me  ' 

and 

C{PMV)  =  {AT  r-'A)-’  .  (31) 

When  Z”1  is  diagonal  and  is  used  as  the  weighting  matrix  W,  the  WLS 
estimate  and  covariance  matrix  in  TRACE  is  also  MV  or  Markov. 
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Nothing  so  far  lias  been  assumed  about  the  actual  form  of  the  distribution  of 
the  random  errors.  If  a  specific  probability  density  function  is  assumed, 
then  it  is  possible  to  seek  the  estimate  that  maximizes  the  probability  or 
likelihood  of  the  resulting  residuals.  In  the  case  of  a  joint  normal  {or 
gaussian)  distribution  of  observational  errors  with  covariance  matrix  Z, 
the  maximum  likelihood  (ML)  estimate  reduces  to  MV.  (In  practice,  as 
Magness  and  McGuire  note  in  Reference  4,  the  gaussian  assumption  is  always 
made.  Reference  4  also  contains  an  excellent  comparison  of  LS  and  MV 
e stimation.  ) 

In  summary,  for  uncorrelated  observational  errors,  the  WLS  estimate  in 
TRACE  is  also  MV;  if  the  errors  are  further  assumed  to  be  normally  distri¬ 
buted,  the  estimate  is  also  ML. 

2.  5.  3  Q- Parameters 

Observations  are  functions  of  many  parameters,  including  six  orbital 
parameters,  differential  equation  parameters  such  as  drag  and  spherical 
harmonic  coefficients,  radar  station  locations,  and  observational  biases.  In 
principle,  all  such  parameters  can  be  estimated,  given  sufficient  observations, 
and  the  covariance  matrix  reports  the  accuracy  with  which  they  have  been  or 
could  be  determined.  In  practice,  however,  only  a  selected  set  of  these  are 
estimated.  There  arises  then,  both  in  actual  orbit  determination  and  in 
systems  analysis,  the  question  of  the  effect  on  (a)  the  parameters  P  being 
estimated,  or  (b)  the  trajectory  itself,  of  errors  or  uncertainties  in  the 
remaining  parameters  O.  The  treatment  here  follows  that  of  Magness  and 
McGuire  in  Reference  5. 

In  TRACE,  the  computation  of  Q-parameter  effects  is  restricted  to  the  error 
analysis  link  FEIGN,  in  which  P-parameter  and  trajectory  covariance  matrices, 
with  and  without  Q-parameter  uncertainties,  are  computed  and  printed.  The 
transformation  from  orbital  (initial  condition)  parameters  to  trajectory 
coordinates  is  covered  in  Section  2.  5.4. 
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Now  the  true  observation  vector  O^JP^,  Q^)  is  a  function  of  true  but  unknown 
values  of  both  P  and  O  and  the  measured  vector  is 


O 

m 


=  °c<PfQt> 


£ 


(32) 


The  vector  of  residuals,  in  the  "measured-minus-computed'’  sense,  is 


O  =  O  -  O  (P  ,Q  ) 
me  m  coo 


°c<PfQt>+e 


-  O  (P  , 
c  o 


Qo> 


(33) 


As  in  Section  2.  5.  1,  we  assume  a  model  for  the  true  observations  O  (P  ,Q  ), 

c  t  t 

which  is  linear  in  a  correction  AP^  to  an  approximate  value  Pq 


°c(PfQt)=  °c(Po-Qt>  +  VPt  •  (34) 

But  now  in  forming  the  "computed"  quantities  we  are  uncertain  as  to  the  true 
value  of  the  Q  parameters  and  must  use  an  approximate  value  Qq  in  the 
calculations.  Thus  O^JP^.Q^)  is  related  to  0^(Pq,  Q^)  by 

O  (P  ,Q  )  =  O  (P  ,  Q  )  -  A  AQ  .  (35) 

c  o  o  c  o  t  q  t 

Collecting  these  results  we  have  the  following  representation  of  the  residual 
vector  O 

me 


O  =  A  AP  +  A  AQ,  +  £  .  (36) 

me  p  t  q  t 

(This  equation  could  have  been  presented  more  briefly;  the  longer  presentation 
is  used  to  make  clear  the  sources  of  the  terms  in  the  residual  vector.  ) 
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t 


I 


) 


Summarizing  the  notation  above, 
observations. 


O  and  O  are  the  measured  and  computed 
me 


P  and  0  =  estimates  of  the  true  parameter  values  P  and  Q 

o  o  r  t  t 

AP.  =  P  -  P  and  AQ  -  Q  -  Q 
t  t  o  t  t  o 

ao  ao 

0  C 

A  and  A  =  the  matrices  and  ,  respectively,  and 

p  q  o\-i.  J 

r  o  o 

€  =  the  vector  of  observational  errors  with  covariance 
matrix  E(t  =  Z. 


The  WLS  problem  is  still  that  of  minimizing  f(AP)  =  |[W*^(0  -A  •  AP)  !  [  ^ 

me  p 

and  the  solution  as  before,  is, 


AP"  =  AP'  =  (aTWA  )_lArW0  (37) 

\  P  pl  p  me  '  ' 

However,  the  covariance  matrix,  as  would  be  expected,  is  affected  by  the 
O-parameter  uncertainties.  Thus 

AP"  =  (aTWA  V*ATW(A  AP  +  A  AQ  +  €) 

\  p  pi  p  P  t  q  t 


or 


oP"  =  AP"  -  AP  =  (aTWA  \“1ATW(A  AQ^  +  0  (38) 

t  \  p  pf  P  q  t 

It  is  seen  immediately  that  if  E(AQ^)  =  0  (meaning  that  unbiased  estimates 
of  the  Q  parameters  are  being  used)  and  E(«)  =  0,  then  E(6P"  )  =  0  and 

AP"  is  an  unbiased  estimate  of  AP  .  The  covariance  matrix, C(P"  )  = 

T  1  - 1 

E(6P"  6P"  )t  is,  by  forming  the  indicated  product,  choosing  W  =  2 

(requiring  uncorrelated  observational  errors)  and  taking  the  expected  value, 

C(P")  =  (aTWA  V  1  +  (aTWA  V>ATWA  C(Q)ATWA  (aTWA  V1  (39) 

\  p  pl  \  p  p/  p  q  q  p\  p  pj  v  ' 
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where  C(Q)  =  E(AQtAQt  )  is  the  covariance  matrix  of  the  Q  parameters. 

From  the  formula,  we  see  that  the  uncertainty  in  the  estimate  P"  =  +  AP" 

as  represented  by  the  covariance  matrix  C(P"  ),  is  in  two  parts.  The  first 

T  -  1 

term  is  CIP'  )  =  (A  WA  }  ,  the  covariance  matrix  (or  uncertainty)  in  the 

P  P 

estimate  P"  resulting  from  random  noise  in  the  observations.  The  second 
term  contains  the  additional  uncertainties  in  the  estimate  of  the  P  parameters 
for  having  used  uncertain  values  of  the  Q  parameters  in  the  process. 
Obviously  C(P"  )  reduces  to  C(P'  )  for  C(Q)  =  0. 

The  effect  upon  the  estimate  P"  =  Pq  +  AP"  of  an  error  (as  opposed  to  an 
uncertainty)  in  a  Q  parameter  can  also  be  predicted.  The  estimate  P"  (Qq) 
of  P  using  the  value  Oq  is 


p"(Q  )  =  p  +  ap"(q  )  =  p  +  (atwa  |'’atwo  (Q  ) 

o  o  '  o  o  \  p  p/  p  me  o 


and  similarly,  using  an  "erroneous"  or  alternate  value  Q 


P"  (Q 


P  +AP"(Q, 
o  1 


=  p  +  (at\va  \'1atwo  (Q,) 

o\p  p  /  p  mcl 

/  t  .  ,  T  f  90 

=  P  +  (A  WA  )‘  A  W  O  (Q  )  +  ■  -JP  C(Q  -  Q  ) 
o  \  p  p  /  p  me'  o'  3Q  1  o' 


Then  the  difference  in  the  estimates  is 


Q  O 

P"(Q,)  -  P"<Qo>  [(ApWAp)‘lApW^fj(°l  -  Qo> 
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a°c 

Since  O  enters  negatively  into  O  (O  =  O  -  O  )  and  since  -  A  . 

c  **  7  me  me  m  c  3Q  q 

we  have 

P"(0.)  -  P"  (Q  )  =  F-(aTWA  V1ATWA  1(0,  -  Q  )  (43) 

l  o  L  \  p  p/  p  qj  i  o' 

The  content  of  the  square  brackets  is  evidently  the  matrix  of  partial  derivatives 
9P" 

aO"  * 

o 

Now  C(P"  )  can  be  rewritten  as 

C(P")  =  c<p'> +  (fer)c<Q>(fcr  ^  <44) 

2.  5.  4  Parameter  Transformations 

The  covariance  matrices  C(P'  )  or  C(P"  )  would  normally  include,  roughly 
speaking,  the  uncertainties  in  the  orbital  parameters  due  to  observational 
errors  and  Q-parameter  uncertainties.  The  orbital  parameters  might  be 
spherical  coordinates  at  time  tQ,  for  example.  Quite  evidently  their  un¬ 
certainties  are  not  very  descriptive  of  the  resulting  trajectory,  period, 
observational,  or  other  related  uncertainties;  hence  the  need  for  transforming 
the  basic  covariance  matrices  to  other  coordinate  systems  and  reference 
times.  A  common  requirement,  for  example,  is  the  covariance  matrix  of 
satellite  position  and  velocity  at  time  t,  resolved  into  radial,  in- track  and 
cross-track  components. 

Suppose  that  a  set  of  parameters  X(t),  intentionally  suggesting 

.  T 

X  =  (x,  y,  z,  x,  y,  z)  at  time  t  *  tQ,  is  related  to  our  P  and  O  parameters, 
and  that 

6X=|^6P+|^6Q  (45) 

o  o 
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This  relation  gives  the  first-order  effect  upon  X(  Pq  ,  Qq,  t)  of  variations 

6P  and  6Q  from  the  nominal  values  P  and  O  .  If  P  is  an  unbiased  estimate 

o  o  o 

P"  of  the  true  vector  Pfl  so  that  E(6P"  )  =  0,  then  the  uncertainty  in  X{t),  due 
to  random  errors  and  Q-parameter  uncertainties,  can  be  written 


jv  _  3  X  D//  3  X  j.  ^ 

6X=lp-6p  +3(3-60 

o  o 


9X.(ATWA  V1AT\V{A  6Q  +  0  + 

\  p  p/  p  q 


ax 


ap0VV“p/  1  '  3"Qo 


6Q 


ax  ap"  ax 

ap  aQ  aQ 
o  o  o 


> 


_  (aTWA  V  ‘  A  A  W c 
3P  \  p  pi  p 


(46) 


Since  we  are  assuming  that  E(c)  =  E(6Q)  =  0  (the  latter  meaning  that  Qq  is  an 

unbiased  estimate  of  Q),  and  that  €  and  6Q  are  uncorrelated  random  variables 

T 

we  see  that  E(6X)  =  0  and  that  the  covariance  matrix  C(X(t))  is  E(6X  6X  )  or 


C(x(t))  ;(^)C(P'  ’fej1 


/ax  ax  ap 
laQ  "  ap  H TT 

\  O  O' 


ax 

aQ 


ax  ap"^ 
ap  KT/ 

n  n  ' 


(47) 


This  is  the  general  formula  for  transforming  P-  and  Q-parameter  covariance 

matrices  into  the  covariance  matrix  for  any  set  of  related  parameters.  Using 

•  T 

the  suggested  interpretation  X  =  (x,  y,  z,  x,  y,  z)  ,  the  partial  derivates 

ax  a  x 

— and  arc  either  just  the  solutions  of  the  variational  equations  (insofar 

as  P(  and  0_^  represent  trajectory  parameters),  or  are  computed  from 
geometrical  relations  (in  the  case  of  obse rvational  parameters).  For  other 
sets  of  parameters  at  time  t,  such  as  tne  spherical  coordinates  R(t)  - 


3R 


3R 


(a,  6,  p,  A,  r,  v)  ,  the  partial  derivative  matrices  v.p-  and  can  be 

,  aR  aR  ax  ......  f  aR  d  °  . 

computed  as  a p  =  -^p—  (and  similarly  for  — ),  wherein 


3P 


3Q 
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3R  2  2 

is  computed  simply  from  the  equations,  such  as  r  =  (x  +  y  + 

which  relate  the  spherical  and  cartesian  coordinates  at  any  time. 

the  covariance  matrix  is  most  easily  obtained  from  C(X(t)}  by 


2,1/2 
z  ) 

In  practice, 


C(R(t))  =  (fx)c(x(tn(fx)T 


(48) 


SECTION  3 


EQUATIONS  AND  METHODS 

3.  1  COORDINATE  SYSTEMS 

There  are  two  coordinate  systems  employed  in  TRACE.  The  earth- center ed 
inertial  system,  known  as  the  "mean  equator  and  equinox  of  date,"  is  basic  to 
all  the  computations,  and  position  and  velocity  in  this  system  may  be  expressed 
in  one  of  three  types  of  coordinates  (paragraph  3.  1.1).  A  station -de pendent 
system  has  al&o  been  introduced  to  facilitate  computations  involving  radar 
observations  and  data  scudies  (paragraph  3.  1.  Z). 

3.  1.  1  Earth-Centered  Inertial  System 

The  basic  coordinate  system  is  as  follows: 


Z 


Figure  Z.  Earth- Centered  Coordinate  System 

where 

0  is  the  center  of  the  earth 

X  is  a  vector  from  0  in  the  equatorial  plane  directed  to  the 
vernal  equinox  at  tCT,  0  hour  GMT  of  launch  date 

Y  is  a  vector  from  0  perpendicular  to  X  in  such  a  direction  that 
(X,  Y,  Z)  is  a  right-handed  system 

Z  is  a  vector  from  0  perpendicular  to  the  equatorial  plane  and 
directed  north. 


The  position  and  velocity  of  a  body  at  a  point  P  may  be  expressed  in  rectangula 
or  spherical  coordinates,  or  in  terms  of  the  classical  (elliptic)  elements  of  its 
orbit,  as  shown  in  the  following  three  paragraphs,  respectively. 

3.1.  1 .  1  Rectangular  Coordinates 

P  =  P  (x,  y,  z,  x,  y,  z)  where  x,  y,  z  are  the  components  of  position  of  the 
body  in  the  X,  Y,  Z  directions,  respectively,  and  x,  y,  z  are  the  components 
of  its  velocity  in  these  directions. 

3.1.  1.2  Spherical  Coordinates 

F  =  P  (a ,  6,  (3,  A,  r ,  v) 


N 


whe  re 


Figure  3.  Spherical  Coordinates 


V  =  a  vector  equal  in  magnitude  and  direction  to  the  velocity  of 
the  vehicle  at  P 


a  =  right  ascension  measured  from  the  X-axis,  positive  eastward 
6  =  geocentric  latitude 

S  =  angle  between  V  and  the  geocentric  vertical  at  P 


a  =  semi-major  axis 


e  =  eccentricity  =  va^  -  b^/a  (b  -  semi-minor  axis) 
i  =  inclination  of  the  orbit  plane 
U  -  right  ascension  of  the  ascending  node 

Li  -  angle  between  the  direction  of  perigee  and  the  line  of  nodes 
t  =  time  in  minutes  from  t  of  last  perigee  passage. 

3.  1.2  Str^ion  -  Dependent  System 
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3 .  2  INITIAL  CONDITIONS 

The  parameters  of  the  orbit  may  be  input  in  any  of  the  coordinates  described 

in  paragraph  3.  1.  1.  The  trajectory  computations  require  earth- centered 

inertial  coordinates;  the  output  includes  spherical  coordinates  and  elements. 

The  formulae  for  the  necessary  transformations  follow.  The  date  chosen  to 

determine  the  X-axis  is  t  ,  zero  hour  GMT  of  launch  date;  the  time  t  at 

g  o 

which  the  parameters  are  specified  is  with  reference  to  this  date. 

3.2.1  Spherical  to  Rectangular 

x  =  r  cos  6  cos  a 
y  =  r  cos  6  sin  a 
z  =  r  sin  6 

x  =  v  [cos  o  (-  cos  A  sin  6  sin  6  +  cos  B  cos  6)  -  sin  A  sin  6  sin  ot~\ 
y  =  v  [sin  o'  {-cos  A  sin  B  sin  6  +  cos  B  cos  6)  +  sin  A  sin  B  cos  o'] 
z  =  v  [cos  A  cos  6  sin  8  +  cos  &  sin  6  1 

If  longitude  (/)  is  input  instead  of  a,  a  is  computed  as  in  paragraph  3.  1.2. 

In  this  case  t  -  t  =  t  . 

g  o 

3.2.2  Rectangular  to  Spherical 


o'  =  tan  ^(y/x) 

6  =  tan-1  ^z  /  V  *2  !-  ^ 

8  =  cos'  1  [  (xx  +  yy  +  zz  )/rv] 


A  =  tan 


1  [~  r(xy  -  yx) _ 1 

[y(yz  -  zy)  -  x(zx  -  xz)J 


r  =  ^  x2  +  y2  + "  z? 
v  =  \]  xZ  +  y2  1  z2 
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3.  2.  3 


whe  re 


Elements  to  Rectangular 


x  =  x  P  +  y  Q 

D)  X  ;0)  X 

y  =■  x  P  +  y  Q 
yj  y  (jj  y 

z  =  x  P  +  y  Q 

u:  z  '  w  z 

x  =  x  P  f*  y  Q  ,  etc. 
<J)  x  o>  x 


P  =  cos  P  cos  u:  -  sin  0  sin  uj  cos  i 
x 

P  =  sin  0  cos  u:  +  cos  0  sin  ujcos  i 

y 

P  =  sin  uj  sin  i 
z 

Q  =  -cos  Qsin  oj  -  sin  P  cos  o'  cos  i 
x 

O  =  -sin  0  sin  uj  -f  cos  P  cos  0)  cos  i 

y 

O  =  cos  Oi  sin  i 
z 

p  =  a (1 -e^),  (semi-latus  rectum) 
u  =  gravitational  constant 

n  =■  =  mean  motion 

M  =  n(t  -  t)  =  mean  anomaly 

E  =  solution  of  (M  =  E  -  e  sin  E)  =  eccentric  anomaly 


r  =  a(l  -  e  co s  E) 
x  =  a(cos  E  -  e) 


0) 


y0,  =  V  I  apl  sin  E 

x  =  -  sin  E 


,-r  _  V  Lip 


0) 


cos  E 


U) 
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(These  formulae  are  for  the  ellipse;  if  the  conic  is  a  hyperbola  (e  >1),  E  is 
the  solution  of  M  =  e  sinh  E  -  E;  and  sin  E  and  cos  E  above  are  replaced 
by  sinh  E  and  cosh  E.  ) 


3.2.4  Rectangular  to  Elements 


- 


e  cosh  E)  -  (e  sinh  E)' 


for  hyperbolic  orbits 


> 
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where: 


r=  V= 2 " ” 


rx  +  y  +  z 


2  .2  .2  .2 
v  =  x  +  y  f  z 


e  cos  E  =  1 - 

a 


„  xx  +  yy  f  zz 
e  sin  E  ~  - — - 

■y/l  a  I  u 


2  2,.  .  .  ,2 

r  v  (xx  f  yy  f  zz ) 

p  _  _ 


e  cosh  E  = 


e  sinh  E  = 


for  the  hyperbola 


D  = 


xx  +  yy  f  zz 


eg. 


D  = 


e  cos  E 
e  r 


H  = 


H  = 


1 


(r  -  P) 


yjlip 

1  xx  +■  yy  +  zz 


yfWp 


P  =  Dx  -  Dx 

x 


Py  =  Dy  -  Dy 

P  =  Dz  -  Dz 
z 

O  =  Hx  -  Hx,  etc. 
x 
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M  =  E  -  e  sin  E  or  e  sinh  E  -  E 

T-,  A  -i/esinE  \ 

E  -  tan  (  - I 

\  e  cos  E  I 
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RADAR  DATA 


3.  3 

All  radar  observation  input  is  converted  to  the  basic  set?  R,  A,  E,  R,  P, 
Q,  P,  Q.  This  section  contains  the  formulae  for  these  conversions.  Input 
is  generally  in  feet,  degrees,  and  seconds,  while  internal  units  are  earth- 
radii,  radians,  and  minutes.  The  details  of  units  conversion  have  been 
omitted. 

First,  however,  Figures  6  through  9  depict  most  of  the  radar  quantities  in 
Sections  3.  3  and  3.4. 


w3  --z 


Figure  6.  Radar  Station  Coordinates 
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Figure  7.  Azimuth  and  Elevation  in  Station  Coordinate  System 

W  is  position  of  the  vehicle 
s 

W"  is  position  of  the  station 

s  s 

Q  is  the  projection  of  Q  =  W  -  W  onto  the  tangent  plane  at  W 


w 


HERE  ARE  PARALLEL  TO  W(  ,  W?  ,  WJt  RESPECTIVELY 

Qpp!S  PROJECTION  OF  0  ONTO  £  .  rj  PLANE 
D;  TOPOCENTRIC  OECLINATION 
HA=  TOPOCENTRIC  HOUR  ANGLE 


Figure  8.  Hour  Angle  and  Declination  in  Station 
Coordinate  System 
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Figure  9.  Station  Network  for  Interferometer  Data 

Here  S,  S^,  and  Sq  are  a  network  of  stations  that  report  range  and  range 
rate  differences.  Let 


R=  |W-S|  ,  Rn=  |W-SJ.  Rq=  |  W-S0  | 


o 


Then 


P  =  R  -  R. 


Q  =  R  -  R 


Q 


P  =  R  -  R, 


3.  3.  I 


Q  =  R  -  R 


Q 


Hour  Angie,  Declination  (HA,  D)  to  A,  E 
_  1 

E  =  sin  (sin  $  sin  D  +  cos  §  cos  D  cos  HA) 


A  =  tan 


■(- 


sin  HA  cos  D 


s  $  sin  D  -  sin  $  cos  HA  cos  D 


where  <P  =  geodetic  latitude  of  the  station. 


3.  3.  2  Right  Ascension,  Declination  (P  A,  D)  to  A,  E 

HA  =  a  -  RA;  compute  a  as  in  paragraph  3.  1.2,  and  then  apply  paragraph  3.  3.  1. 
3.3.3  Af,  At,  to  R,  R 

R  =  M  where  is  input 

R  =  k^  At  where  k^  is  input 
3.  3.4  L,  M,  N  to  A,  E 
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PARTIAL  DERIVATIVES  OF  RADAR  DATA 


3.  4 

In  tracking  and  data  studies,  it  is  necessary  to  compute  partial  derivatives  of 
radar  data  with  respect  to  parameters  of  the  initial  conditions,  differential 
equations,  station  locations,  and  observations.  For  the  purposes  of  this 
section,  it  will  be  assumed  that  the  (integrated)  position  of  the  vehicle  is 
known,  in  earth  -  center  ed  inertial  rectangular  coordinates,  as  are  the  partial 
derivatives  of  these  coordinates  with  respect  to  the  first  two  types  of 
parameters. 

3.4.1  Notation 


p. ,  i=l,  2,  ...,  n 
ri 


bx  by  bz 

’  ^?i  ’  '  "  ’  ^p. 


I 


h 


the  ordered  list  of  initial  condition  and 
differential  equation  parameters  for 
which  partial^  are  to  be  computed 

the  partial  derivatives  of  x,  y,  .  .  .  ,  z 
with  respect  to  the  pi 

longitude  of  the  station 

geodetic  latitude  of  the  station 

height  of  the  station 


O' 


w. , 
J 


w. 

J 


1  =  1. 


2,  3 


s  s 
Wl*  W3 

e 


a 

e 


b 

e 


ae(l-s) 


o  +  u;  (t  -  t  )  +  as  in 
g  e  g 

position  and  velocity  of  the  vehicle  in 
the  station-dependent  W-system 

position  of  the  station  in  the  above  system 

ellipticity  of  the  reference  ellipsoid 

semi-major  axis  of  the  earth 

semi-minor  axis  of  the  earth. 


I 
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3.  4.  2  Position  and  Velocity  in  the  W-system,  and  Associated 

Preliminary  Computations 

Wj  =  x  cos  a-  -I-  y  sin  ry 

w ^  =  -x  sin  <y  -I-  y  cos  <-k 

=  7. 

Wj  =  (x  f  '}'ey)  cos  Of  +  (y  -  'J>ex)  sin  o' 

=  -(x  +  m:  y)  sin  a  +  (y  -  u;  x)  cos  a 
e  e 

=  z 

To  transform  the  partial  derivatives  of  the  Earth- Center ed  Inertial  (ECI) 
rectangular  coordinates  to  the  station-dependent  system,  the  following 
quantities  are  necessary. 

Differentiation  of  the  above  six  equations  shows  that  a  simple  substitution  of 


Sw.  3w.  > 

fo r  w . ,  for  w.,  j  =  1 ,  2 .  3 ,  and  —  ,  . 

5p-  J  3Pj  J  ^Pi 


‘  *  ’  *P; 


for  x,  y . z  yields 


t)w 

1  ay 

-r -  =  —  cos  o  +  ^r~  sin  rv 

5p.  >>p.  5p. 


ow 

2  -ox  .  ,  dy 

-  =  - - sin  of  cos  O', 

'‘n  An  -'n 


ip.  bp. 
1  1 


Differentiating  with  respect  to  Si,  since  =  1,  gives: 


=  -X  Sin  Of  +  y  COS  ry  =  w. 
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— —  -  -xcos  o  -  y  sin  a  =  -w^ 


bT~  ~  w2 


-  =  -  w , 

1 


To  find  the  station  position  in  the  W-system,  we  use: 


where: 


w,  =  (a  A  +  h)  cos  $ 
1  e  s 


w^  =  (b  B  +  h)  sin  $ 
3  '  e  s 


A  =  (cos2  $*  +-|  sin2  $:V1/2 
s  ^ 

ae 


B  =  (sin2  i*  +  -§•  cosbV1'2 

S  b2 

e 


Differentiating  with  respect  to  and  h 


5wl  s  <’> 

^  =  ‘W3 


s  (1) 

5"  =  w! 

3$  1 


Approximate  formulas;  correct  only  for  a  spherical  earth. 
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At  this 
and  V, 


$w' 


=  cos 


dh 


point  it  is  convenient  to  introduce  three  intermediate  vectors  Q,  U, 
and  the  quantity  R^. 

g 

O  =  W  -  W‘  vehicle  position  relative  to  the  station. 
q}  =  Wj  -  w* 

q2  =  w2 
q3  ,  w3  -  w* 

R  =  i  Q|  =  ^qf+  qf  +  q3 

U  =  Q/R,  a  unit  vector  in  the  direction  of  Q 
Uj  =  qx/R 
U2  =  q2/R 

u3  =  q3/R 

V  is  the  vector  U  referred  to  the  East-North-Up  system 
V1  =  U2 

v^  =  -Uj  sin  $  +  u^  cos  § 
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cos  I  +  sin  $ 

nr  t 

v=  Vvl  +  V2 


Then, 


Rj  =  vR 

=  sin  E 


v  =  cos  E 


v 


=  cos  A 


—  =  sin  A 
v 


To  compute  the  and  we  will  need  v-,  and  to  this  end  we  compute 

u  =  I  [w  -  ur3  =  [w  -  (u  .  w)  u] 


then 


V\  =  U2 


=  -u^  sin  $  +  u^  cos  $ 


v  = 


V1V1  + V2VZ 


3.  4.  3 


Range  Partials 

V~~2  2  2 

q j  +  +  ^3  results  in: 


SR  „  30  *ql  ,  "*Jq2  .  3q3 

dp.  dp.  1  dp.  2  dp.  3  dp. 
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3.4.  5 


dA  1  .  .  . 

ar  =  -r  (v2vi  -  vivz) 

v 


Elevation  Partials 


Differentiating  E  =  sin  ^  v^  =  cos  ^  v, 


:  j_  (  awi 

rRi  \aPi 

(' 


*  ,  dw3  .  *  oR 

cos  $  +-  -  sin  $  -  t —  sin 

§p.  ap7 


1  f  **  ■  **  ^R  •  ~ 

w,  cos  $  -  w  sin  $ - —  sin  E 

a*’  *vi  V  a#' 


BE  1  .  *  SR  .  . 

SI“S*7  w2  COS  *  '  9T  SinE) 

BE  _  -1  ..  ,  dR  . 

dh”  =  R^  (1+3h  smE) 


be 

3E 


bias 


1 


BE 

}t 


u  ^  cos  $  *  +  sin  $  " 
cos  E 


3.4.6  Range  Rate  Partials 

Differentiating  R  =  (U  •  W  ), 


) 
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ff  =  «w2di  -  wlu2>  +  (^2U1  -  Wlu2> 


aR 

ar  =  -  ui  cos 


-  sin  § 


aR 

aR.. 

bias 


=  1 


aR 

at 


=  R  =  {U  *  W)  +  (U  ■  W) 


where  W  =  - 


(“e  j  + 


W  +  2  LX 


J  = 


-x  x  sin  (y  +  ij‘  y  cos  a 
/  e  e  ’ 


LX  = 


-x  x  cos  O'  -  o)  y  sin  or 
e  e 


( 


1/2 


3.4.7  P,  Q,  P,  Q,  Partials 

These  partial  derivatives  are  obtained  by  differencing  the  R,  R  partials  using 
the  appropriate  station  locations. 
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3.4.  8 


A,  E  Partials 
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3.  5.  3  Observations  with  Normally  Distributed  Random  Noise 

o  =  o  4-  r  o  is  the  noisy  observation  {of  type  j  from  station  s) 

c  n 

o^  is  the  nominal  computed  observation 

r  is  the  noise  added, 
n 

r  =no  +6  a.  is  the  appropriate  sigma  for  type  j,  station  s 
n  sj  sj  SJ 

3  .  is  the  appropriate  bias  (if  any) 
sj 

n  is  a  random  element  from  a  set  of  numbers 
with  mean  zero  and  unity  standard  deviation. 


I 

9 


3-25 


TRAJECTORY 


The  position  and  velocity  components,  X  =  (x,  y,  z)  and  X  =  (x,  y,  z),  of  the 
vehicle  and  their  partial  derivatives,  XD  and  X  ^  (i  =  1,  .  .  .  ,  n),  with 
respect  to  the  trajectory  (initial  condition  and  differential  equation)  param¬ 
eters  are  functions  of  time  defined  by  their  differential  equations  and 
appropriate  initial  conditions  (paragraphs  3.  6.  1  and  3.  6.  2).  The  equations 
are  integrated  numerically  (paragraph  3.  6.  3),  and  at  each  observation  or 
print  time  all  the  quantities,  X,  X,  X^.,  and  X^  (i  =  !  ,  .  .  n),  are 

obtained  by  interpolation  (paragraph  3.  6.4)  in  the  integrated  results;  from 
these  the  computed  radar  observations  and  partial  derivatives  (Sections 
3.  3  and  3.  4)  and  the  trajectory  output  (paragraph  3.  6.  5)  are  computed. 


3.  6.  1 


Differential  Equations 


The  equations  of  motion  of  the  vehicle  are 


X  =  ^*  +  F 


where  p  is  the  gravitational  constant  (GM)  of  the  earth,  r  ]x(  = 

(xZ  +  yZ  +  and  F  =  F^  +  F  ,  +  F^  is  the  perturbative  acceleration 

dje  to  asphericity  of  the  earth,  extra-terrestrial  gravitational  forces,  and 
atmospheric  drag,  respectively.  The  initial  conditions  X(t  )  and  X(t  ), 
if  not  given  directly,  are  computed  from  the  initial  spherical  coordinates 
or  elliptic  elements.  See  Section  3.  2  for  these  formulae. 

The  perturbative  acceleration  F(  due  to  the  asphericity  of  the  earth  is  derived 
from  the  assumed  potential  function. 


5 

.  Ji  4  n  /a 

xn 

U  rH 
r 

■  Z’i 

—  )  Pni(  sin  6)  cos  m(X  -  X  ) 

/  n  nm 

n=  2 

n = 2  m  =  1 

_ 
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The  force  vector  in  the  ECI  coordinate  system  is  then: 


cos  cp 

COS  O' 

-  sin  0 

-sin  cp  sin  0/ 

1  BU 

cos  cp 

sin  0 

cos  a 

-sin  cp  ain  <y  j 

[  bE 

sin 

0 

cos  cp  ' 

^  % 

where  a  =  0/  +  <y  (t  -  t  )  is  the  right  ascension, 

g  e  g  B 

The  gravitational  attraction  of  other  bodies  contributes 

k  /  y  v  v  v 


K  /  X  -  X.  X.  \ 


where  m.  is  the  mass,  relative  to  the  earth,  of  the  jth  body  and  X.  is  the 

J  th  ^ 

vector  position  of  the  j  body,  as  obtained  from  the  JPE-STE  planetary 

coordinate  tapes.  For  a  description  of  these  tapes  and  their  preparation  see 

Reference  6. 

Note  that  the  tabular  planetary  coordinates  are  with  respect  to  the  Mean 
Equator  and  Equinox  1950.  0  coordinate  system,  whereas  TRACE  calculations 
are  referred  to  0  hour  GMT  of  start  day.  The  planetary  coordinates  are 
transformed  to  the  coordinate  system  of  TRACE  before  is  calculated. 

The  subroutine  is  described  in  Reference  7  which  in  turn  refers  to  Reference  8. 

The  effect  of  atmospheric  drag  is  the  term 


_  va/cda\x 

F3  =  PtHHXA 
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where  o  is  the  density  at  height 


h  =  r  - 


aeU-€) 


2  2 

,  2>  x  +  y 

1  -  (2  0  -  e  )  - ^ — 
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above  the  oblate  earth,  and  where 


cda 

w 

is  the  dr. 

*A 

is  the  ve 

atmos  ph( 

*A  = 

x  +  j:  y, 
e7 

ii 

< 

V  -  y  x, 

7  e 

±A  = 

z  ,  and 

< 

> 

ii 

lxA|. 

is  the  drag  coefficient,  (or  "ballistic  coefficient"), 


The  atmospheric  density  is  computed  from  an  atmosphere  model  (or  certain 
combinations  of  models)  given  by  References  9,  10,  and  11. 

3.  6.  2  Trajectory  Partial  Derivatives 

The  partial  derivatives  of  vehicle  position  and  velocity  with  respect  to 
trajectory  parameters  can  be  approximated  analytically,  or  can  be  obtained 
by  a  simultaneous  numerical  integration  of  the  variational  equations. 

3.  6.  2.  1  Variational  Equations 

The  variational  equation  for  an  initial  condition  parameter  a  is 


X 

Of 


with  initial  conditions  X  (t  ) 

y  o 


X  (t 

O'  o 


(54) 
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For  a  differential  equation  parameter  8. 


X„  = 


‘  a 

/  UX  N 

,  3f1 

ax 

)  + 

x*  + 


5F 

ax 


aF 

ap 


55' 


with  Xp(to)  =  Xp(to)  =  0. 

Here  X  =  -l~ »  X  =  -§~~»  X_ ,  X  and  are  all  3-vectors.  The  contents 
a  Da  cl  9a  (3  (3  9(3 

3  F 

of  the  square  brackets  and  are  3x3  matrices.  The  system  is  solved 
for  each  parameter,  and  all  the  numerical  integrations  are  carried  out 
simultaneously. 


9F 

In  the  above  equations  the  principal  contributions  to  ^  stem  from  the 

oblateness  coefficient  and  from  the  dependence  of  the  drag  force  upon 

the  position  of  the  vehicle.  {The  latter  is  important  for  low-altitude  satellites.  ) 

Lesser  sources  are  the  other-body  gravitational  forces  and  the  higher  order 

3  F 

oblateness  terms;  they  are  ignored  in  the  calculation  of  — ^  . 


The  matrix  in  square  brackets  is  calculated  as  the  sum  V  +  T  where  V 
derives  from  the  gravitational  force  including  the  oblateness  term,  and 
T  from  the  drag  force.  The  spherical  earth  contribution  is  easily  derived: 


a 

ax 


-u. 


-  3  ax  0  -4  ar 
r  AX  -  3r  x  ax 


(56) 


since  =  — — .  The  oblateness  component  is  not  so  simply  obtained.  It 

is,  of  course,  derived  from  the  contribution  of  the  term  to  the  perturba¬ 
tive  acceleration  Fj,  which  is  in  turn  the  gradient  of  the  potential  U.  The 
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calculation  is  tedious  and  only  the  final  matrix  V  is  given  here. 


P  -OfJ0a2(x2S-U) 

P  +xyJ„a2S 

P  f xz  J_a2T 

xx  2  e 

xy  7  2  e 

xz  2  e 

P  f  xy  J?a2S 
xy  2  e 

P  -Of  J0a2(y2S-U) 
yy  2  e  7  ' 

P  fyz  J_a2T 
yz  7  2  e 

P  Ixz  J ~a2T 

P  fyz J0a2T 

P  -  Of  J0a2(z 

xz  2  e 

y  z  7  2  e 

zz  2  e 

where  is  the  principal  oblateness  coefficient, 


xy 


=  Q  = 

c.  .  w 


3  r 


s  , -Xj  (i 

2r  r 


A  <3  -  > 

2r  r 


2r  r 


The  T  matrix,  which  shows  the  dependence  of  the  drag  force  upon 
position,  is  derived  as  follows 


T 


(  D 


VAXA> 


VAXA  |x  + 


aVA  „  **a) 

ZX-  +  cVA  ~Tx) 


The  derivatives  of  p, 


and  X^  are: 


/p\  -  5p 

y  '  hX  ~  Sh  SX 


(57) 

T-3U) 


the  vehicle 


(58) 
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where 


h  =  r  - 


ae(1  '  e) 


1  -  (2  e  -  e 


2,  x2  ,-2^2 


dh 

ah 

ah 

ah  \ 

ax 

ax' 

ay 

’  az  ) 

dh  _ 

X 

{■ 

aee  (2-3e+e2)  z2 

Bx 

r 

Tr 

T~Z  277“2  2.-1 3/2 

“(2  c  -  g  )(x  fy  )] 

ah  _ 

y 

{■ 

aee(2-3efe2)  z2 

ay  “ 

r 

’ 

ir~z  2. ,  2,  2. -,3/2 

r  ■  (2  e  -  e  )(x  +y  )  J 

ah 

z 

{■ 

_L  _ 

a^  e  (2-3e  +  e  2)(x2+y2) 

az  " 

r 

1 

[r 

2  ,-y  2.,  2  2. -i3/2 

-  (2  e - e  ) (x  +y  ) J 

! 


and  ,  the  rate  of  changes  of  density  with  altitude,  depends  upon  the  model 
atmosphere,  its  parameters, and  h. 


An  approximation  to 


in  the  form 


9P  _ 
3h  " 


P 


-  P 

K 


may  be  used  by  specifying  a  value  for  p'  in  each  of  the  intervals 
0  <h  <  108  n  mi  and  108  <  h  <  378  n  rni.  Alternatively  (as  in  Reference  12), 

g 

may  be  calculated  from  density  expressions,  for  76  <  h  <  108  n  mi 


pl 


=  5. 606  X  10 


-12 


*  »■ 


4/3 


Fio.  7K1  + 


h  -  76/l  +  cosV\3| 

T5  3"  \  2  )  I 
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and  for  108<h<  378  n  mi, 


p2  =  po(h)  (0.  85F10  7)|l  +  0.19  [exp(0.0102h)  -  1.9]  (*  +  c°s.f.)3 

where  logjQ  P0  (h)  =  d7  -  0.  00368  h  +  6.  363  exp[  -  0.  0048h]  . 
Differentiating  each  of  these  expressions  with  respect  to  h,  one  obtains 


and 


3pZ  _  0.00368  +  0. 0305424  exp[ -0.0048  h]  , 

~Zh  =  ‘  p - 0.l34Z9i?8T9 -  +  (0.  dbJFfo  7 

exp  2.302585  (d2  -  0.  00368h  +  6.363  exp[  -0.0048h])j  X 

jo.  001938  (— V— ^  exp  [0.  0102n]j 

r)V 

<b>  -59C  =  ^  [  (x  +  UJey)2  +  (y  -  o.-ex)2  +  z  1  11 


0) 


» 
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Combining  these  results,  the  elements  T^.  of  the  T  matrix  are 


cda 


V  X  ^  ,uexAyA 

W  l  A  A  }h  dx  V  A 

\  A  j 


12  "  2  W 


„  ■  3o  3h,  p  e  A  , 

VAXA  ah  3^  +  "  -  +  PVA“>< 


C°A  ( v  *  i°  3h 
~TW  l  A  A  dh  dz 


cda  / 


(  VA*A 


dp  dh  P0)e^A 
dh  dx  "  VA 


“  PVao;. 


[  v  v  »>  + 

2  W  \  Ay  A  dh  3y  + 


^0  dh  ,  pw;eXA^A 


(v  y  3“  ^ 

2  W  I  AyA  dh  dz 


cda  /„  .  a»  ah  P-^A^A 

2 W  I  VAZA  dh  dx  '  VA 


.  cda  /  ap  ah  pjle*A*A 

1  32  ~  '  2 W  l  VA  A  dh  dy  VA 


T  =.^(v  z  i£  » 

33  2 W  l  A  A  dh  dz 


dF 

The  matrix  — r  is  simply 

ax 


i  „  cda  (  xaxA  .  , 
—  -ijVa  + 1 


3.  6.  2.  2 


Variational  Equation  Initial  Conditions 


The  initial  conditions,  X  (t  )  and  X  {t  ),  are  given  here  for  three  types  of 

a  o  a  o  b  K 

pa  ramete  r  s. 


For  parameters  in  rectangular  coordinates  of  the  initial  position  and 
velocity, 


I,  the  3x3  identity  matrix  and 


For  the  spherical  coordinate  parameters, 
o  (right  ascension) 


dx 

*nr 


-y 


do 


x 


_ 

^  - 


X 


dz 

do 


0 


*  (declination) 


- — 7-  =  -r  sin  6  cos  o 


Sy  _ 

d* 


-r  sin  b  sin  o 
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3z 

—  =  r  cos  6 
oo 


3x 

36 


-z  cos  cy 


3y  _ 

36  " 


-  z  sin  O' 


3z 

'36 


v  (cos  0  cos  6  -  cos  A  sin  8  sin  6) 


8  (flight  path  angle) 

3x  _  3y  _  _3z  _ 

38  “  3^  ~  30  =  U 

=  -v  [(sin  0  cos  6  4-  cos  A  cos  8  sin  6)  cos  a  +  sin  A  cos  8  sin  ol 
3p 

-  -v  [(sin  8  cos  6  +  cos  A  cos  0  sin  6)  sin  o  -  sin  A  cos  p  cos  o'] 

on 

3z 

•“  =  v  (cos  A  cos  0  cos  6  -  sin  8  sin  6) 
do 

A  (azimuth) 


3x  _  3y  ^  3a  _ 
3A  AS."  3A" 


3x 

3A 


v(sin  A  sin  6  cos  a  -  cos  A  sin  o)  sin  8 


3y 

3  A 


v(sin  A  sin  6  sin  o  4-  cos  A  cos  o)  sin  p 
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=  -v{sin  A  cos  6  sin  9) 


r  {magnitude  of  radial  vector) 


dx  _  x 
dr  _  r 


ly  _  y 

dr  r 


dz  _  z 
dr  =  r 


dx  _  dy  _  dz  _  n 
dr  ~  dr"  “  dr 


v  (velocity) 


dx  _  dy;  _  dz 
dv  “  dv  '  dv 


dx  x 
dv  '  v 


Sy  =  y 

dv  v 


dz  _ 

^v  ~  v 


The  equations  for  the  partial  derivatives  of  position  and  velocity  components, 
with  respect  to  elliptic  elements,  are  used  to  compute  initial  conditions  at 
time  tQ  of  the  variational  equations  for  the  parameters  of  this  type.  They 
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may  also  be  used  to  estimate  analytically  the  trajectory  partial  derivatives. 
These  equations  are  as  follows. 


a  (semi-major  axis) 


X  =  I  (X  -  X  ) 
a  a  2n  ’ 


x  =  I  (X  -  ™  x  -  4  X ) 

a  a  '  2n  2 


X  =  - 


p-X 


e  (eccentricity) 


X  =  -  a  + 


y2 

10 


1  r  ( 1  -  e  )  J 


x  y 

P  +  o 


r(l -e  ) 


X  =  - 
e 


r  (1  -  e 


TTTZ 


y<r 


y  +  n  (— )  x  y 

’  u>  r  <j)  o: 


rd-e2,l/2 


(1-e 


2  _ 

w  .  ,  ,  a  .  2 

2.1/2  y>  r  m 


O 


i  (inclination) 


X.  = 


1  (P2  l-  q2)1/2 

z  z 


W 


wlie  re  W  =  P  x  O 


X.  =  - « - ~-y~TT?  W 

1  (p2  +  oV7^ 
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Q  (longitude  of  ascending  node) 


5x 

50  =  -y 

50 

sy  _  x 

5y 

50 

50 

5z  _ 

5z 

50  =° 

50 

o,  (argument  of  perigee) 


X  =  -y  P  f  x  Q 

ui  u:  y, 

X  =  -v  P  4-  x  Q 

U'  '  u 

t  (time  of  perigee  passage) 


X.  =  -X 

y  J.  X 

XT  ■  ~~T~ 

r 

The  variational  equation  l'or  initial  time  is  like  Eq.  (54)  with  the  initial  conditions 

Xt  (t  )  --  -X(t  )  and  X.  (t  )  =  -X(t  ). 
t  o  o'  t  '  o'  '  o' 

c  ° 

3.  6.  2.  3  Differential  Equation  Parameter  Non-homogeneous  Terms 

3F 

The  non-homogeneous  terms  for  the  differential  equation  parameter 
variational  equations  are: 


cda 

w 


(drag  coefficient) 
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u  (gravitational  constant) 
dF  _  Fi  +  F2  _  X 


J.,  J..  ,  X..  (oblateness  parameters) 

1  lk  ik  ^ 

Denote  the  perturbative  force  components  in  the  Up,  East,  North  system 
(see  paragraph  3.  6.  1)  as  follows: 


»U  -  '  “Z 


5  4  n 

LA  +  X]  B  cos  m  (X-X 

,  n  '  nm  i 


n  =  2 


n  =  2  m  =  l 


>E  "*  2 

r 


4  n 

F  C  sin  m  (X-X 

nm 

n=2  m=l 


nm' 


gN  ~  2 


5 

£  D 


+ 


4 

£ 


n 


n=2  n=2 


ll 

EE  cos  m  (X-X 
.  nm  ' 

m  =  1 


nm 


then 


agrT  a. 

-u  i 


aj. 


7  V  SJi 


=  0, 


dgN  -ll  °i 


aj. 


T  T7 

r  i 


(60) 


agTT  B.,  cos  k  (X-X..  )  bg,-,  C..  sin  k  (X-X  ) 

6U  -  ll  ik  lk  ,  feE  -ll  ik _ ik' 


aj. 


ik  r 


J. 


ik 


aj.,  “  ~T 

ik  r 


J. 


ik 


(61 


a§N  Eik  cos  k  (X-Xik> 

aj.,  2  j. 


ik  r 


ik 


(62) 
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(63) 


~2  kBjk  ^nk  (X-X.k), 


kC..  cos  k  (X-X.,  ) 
2  lk  lk 

r 


N 


Sr 

^1N  -u  _ 

a\T  =  ~z  kEik  sink  (x-xik> 

ik  r 


(64) 


The  component  terms  are  then  rotated  to  the  ECI  system  by  the  matrix  given 
in  paragraph  3.  6.  1. 

3.  6.  3  Integration  Methods 

For  the  numerical  integration  of  the  differential  equations  described  in 
paragraphs  3.  6.  1  and  3.  6.  2,  a  choice  of  methods  is  offered.  They  are  widely 
known  as  the  Adams -Moulton  and  the  Gaus s -Jackson  methods,  and  the  sub¬ 
routine  names  are  AMRK  and  DE6F,  respectively.  Both  are  variable -step- 
predictor-corrector  methods  with  automatic  local  truncation  error  control 
and  double  -  preci  sion  accumulation  features.  Both  use  the  Runge-Kutta 
methods  to  obtain  starting  values.  (See  Reference  13.  ) 

The  Gaus  s -Jackson  method,  utilizing  differences,  is  of  higher  order  and 
has  proved  to  be  remarkably  effective  in  the  integration  of  most  satellite 
trajectorie s.  In  some  restricted  but  well-controlled  tests,  this  method, 
applied  to  the  equations  of  motion,  produced  results  that  compared  favorably 
in  both  speed  and  accuracy  with  more  sophisticated  special  perturbation 
methods.  (See  Reference  14.)  Because  of  its  procedure  for  changing  the 
step  size,  the  subroutine's  efficiency  will  drop  and  lose  accuracy  when  the 
step  size  changes  are  extreme,  as  in  highly  eccentric  orbits  or  upon  entry 
into  the  atmosDhere.  In  these  cases  the  use  of  AMRK,  in  which  the  variable 
step  is  more  stably  handled,  is  recommended. 


3.  6.  4 


Interpolation 


Each  time  the  position  and  velocity  (and  their  partial  derivatives)  of  a  vehicle 

are  required,  the  desired  quantities  X,  X,  X  and  X  are  obtained  by 

^i 
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interpolation  from  the  results  of  the  integration.  This  technique  permits  an 
uninterrupted  numerical  integration,  is  comparatively  rapid,  and,  as  used 
here,  is  quite  accurate.  In  particular,  function  values  and  their  first  and 
second  derivatives  at  the  two  adjacent  integration  steps  are  retained  to  permit 
5^  and  3*^  degree  interpolations  for  position  and  velocity,  respectively. 


The  method  used  is  Hermite  interpolation  (Reference  15). 

3.  6.  5  Trajectory  Output 

The  position  and  velocity  vectors,  X  and  X,  of  the  vehicle  are  the  basis  of 
the  trajectory  output.  From  these  quantities,  obtained  by  interpolation  from 
the  results  of  the  numerical  intergration,  are  computed  the  spherical  coordi¬ 
nates  a,  5,  (3,  A,  r,  v  of  the  vehicle  (see  paragraph  3.  2.  2)  and  also 


geodetic  latitude,  $  =  tan 


-1 


T2  2.  1  /2  7z 

(x  +  y  )  (1  -  e) 


longitude,  l  =  a  -  u  -  o;  (t  -  t  ) 

§  e  8 


height,  h  =  r  - 


aeU  -  e 


1  -  (2e  -el 


2,  x2  +  y2 
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These  results  are  output  in  feet,  degrees,  and  seconds. 

The  partial  derivatives,  X  and  X  (i  =  1  ,  .  .  .  ,  n),  of  vehicle  position 

Hi  *"i 

and  velocity  with  respect  to  the  n  trajectory  parameters  can  also  be  printed. 


Optionally,  the  elements  of  the  osculating  ellipse  are  output.  Included  are 
the  elements  a,  e,  i,  ft,  u>,  t  (see  paragraph  3.  2.  4),  and  also 


Mean  anomaly,  M  =  E  -  e  sin  E,  where  E  =  cos 


(deg) 
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(Formulae  from  Reference  16.  } 

3.  6.  6  Analytic  Trajectory 

On  option,  an  analytic  orbit  can  be  obtained  instead  of  an  integrated  orbit. 
The  analytic  trajectory  consists  of  a  Keplerian  orbit  with  nodal  regression 
and  element  decay  due  to  atmospheric  drag.  The  changes  in  elements  are 
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he  se 


calculated  at  perigee  on  every  n  revolution  where  n  is  specified  T 
formulae  are  taken  from  Reference  17. 

3.  6.  6.  1  First  Order  Nodal  Regression 

For  one  revolution,  for  small  e: 


AO  = 


-3tt  J_a  cos  i 
2  e _ 

2  n  zTz 

a  (1  -  s  ) 


Aw  = 


3-n  J.,a2  (4-5  sin2  i) 
c  e 

“  9  2  /i  2v2 

2a  (1  -  e  ) 


or 


n  T  ,2  1/2 

An  -3~Ve  u 

■ST  "  .  7/2“  272 

2a  (1  -  e  ) 


cos  i 


T  2  1/2 
.  3J  a  \i 
Ay  2  e 

W  ~  “7/2“  2.2 

4a  (1  -  e  ) 


(4  -  5  sin2  i] 


3.  6.  6.  2  Element  Decay  Due  to  a  Rotating  Atmosphere 

Define  scale  height  H  =  h  +  12  where  h  is  the  altitude  in  nautical  miles. 


cl0 

If  ~tt  >  2,  the  following  formulae  are  used: 

il 

2 


Aa  =  -  0 


1  f 


1  -  8e  f  3e‘ 

8c  (1  -  e2)  J 


Ae  =  -  O  ( 


1  -  e 


1  - 


3  4  e  -  3e‘ 

8c  (1  -  e2) 
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n 

—  =  ratio  of  earth  rotation  rate  to  satellite  mean  motion 
n 


pp  =  density  at  perigee 


n  ^DA  ne  ,1/2  n  .-1/2 

B  =  •"  -  —  a  o  f  (Zttc) 

m  n  p 
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If  = 
H 


c  <  2  then  use : 


G  (1  +  e) 

Aa  =  -  ^ - t -r*-  [(1  -  2e)  I  (c)  +  2el  (c)] 

(1  -  e)1'*  °  1 

Ae  =  -  x  (0  -  e)  I^c)  +|[lQ(c)  +  I2(c)]} 

Ai  =  -K  rIQ  (c)  -  I2(c)l  +  cos2  a  [I  (c)  -  2el^(c)]}  sin  i 


AO  =  -  K  fl2(c)  -  2el^(c)]  sin  ir  cos  y : 


Aa  =  -  AC  cos  i 


where 


CnA 

-i  D  2  ,  ,  -c . 

G  =  2-  -  a  o  f  (e  ) 

m  p 


C_A  n 

K  =  tt  -  —  a  o  JT~  {e  C) 

m  n  p  ^ 


1^,  1^,  1^  =  imaginary  Bessel  functions  of  the  first  kind. 


3.  6.  7  Initial  Condition  Derivation  (Gaussian  Method) 

On  option,  initial  conditions  may  be  calculated  from  two  sets  of  RAE 
observations.  The  following  procedure  is  adapted  from  Reference  18. 

Let 

=  cartesian  vector  associated  with  first  RAE  observation 
=  cartesian  vector  associated  with  second  RAE  observation 
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g7  =  g2  "  g4 
g8  =  Vg^ 

g9  =  (g7)2 


1  1 

gy  ' 

'l- 

gi  ^ 
g2  } 

1  +  (g^  _  83^4* 

g3  | 
89  ' 

r*  1 

^  2g7 

-■) 

1  -  --  Cg5  +  3(g(l)-  g3)l 
s6 

then  g^1+  ^  =  g^  -  Ag. 
Iterate  until  |Ag|  <  e. 


The  Keplerian  elements  are  then  given  by  the  following  equations. 


X 


+  |X2|  -  2  ^IXji  |X21  cos  g  cos  f 


,  .  2 
2  sin  g 


IV  V 

e  cos  E,  =  1  -  -  ;  e  cos  E0  =  1  -  - 

la  2  a 


e  sin  E^  =  |g  f  e  cos  E^  -  (cos  2g)(e  cos  E^)] 


-(sin  2g)(e  cos  E^) 


E 


1 


_  i 

=  tan 


e  sin  E^ 
e  cos  E^ 


0  <  E^  <  2n 
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=  [(e  cos  E^)2  +  (e  sin 


T  =  t.  -  (E.  -  e  sin  E.) - 

1  1  ^ 


if  a  >  0 


=  cos  *  W 


0  <  i  <  n 


.  W 

-1  x 


0  <  Q  <  2tt 


\/l-eZ 


cos  E^  -  e  vi_e  sinE 

1  -e  cos  E^  ’  sin  V1  1-e  cos  £~ 


U1  cos  Vj  -  sii. 


U1  sin  +  Vj  cos  v, 

f  -1  Pz 
tan  3- 

z 


0  <  uj  <  2r 
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DIFFERENTIAL  CORRECTION  AND 


3.  7 

ASSOCIATED  COMPUTATIONS 

Basically,  the  problem  of  differential  correction  is  the  change,  or 
correction,  of  a  given  set  of  parameters  so  as  to  achieve  some  specified 
result.  In  this  case,  the  goal  is  to  minimize  the  weighted  sum  of  the 
squares  of  the  differences  between  the  (observed)  radar  data  and  the 
corresponding  quantities  computed  from  the  observational  model.  The 
model,  of  course,  includes  the  trajectory  and  radar  parameters  to  be 
corrected. 


3.  7.  1  Notation  and  Nomenclature 

In  general,  matrices  and  vectors  will  be  denoted  by  Roman  capitals; 
their  components  by  corresponding  lower  case  letters,  with  subscripts 
where  applicable. 


n  number  of  observed  quantities 

m  number  of  parameters 

k  number  of  effective  parameters  (m  minus  the  number  of 
restraint  equations) 


o. 

l 


O 


me 


P 

AP 


it^1  observation  (may  be  range,  azimuth,  elevation,  range 
rate  P,  O,  P,  Q) 

radar  sigma  (multiplicative  weighting  factor)  to  be  applied 
to  data  type  j  from  station  s 

radar  bias  (additive  weighting  factor)  to  be  applied  to  data 
type  j  from  station  s 

vector  of  weighted  residuals  (differences  between  observed 
and  computed  radar  quantities) 
vector  of  parameters 
correction  vector  for  P 
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vector  of  bounds  on  solution  AP 

3o. 

matrix  of  weighted  partial  derivatives;  a..  =  - /  <r 

lJ  Pj 

i  =  1,2,  ....  n;  j  =  1,  2,  .  .  .  ,  m 
o.  determines  the  appropriate  <r  to  be  applied 

A  transpose 


constraint  matrix 


3.  7.  2 


Sigmas  and  Biases 


Usually,  some  empirical  evidence  is  available  as  to  the  relative  accuracy  of 
various  types  of  data  from  different  stations.  For  this  reason,  two  types  of 
factors,  each  correlated  with  a  particular  station  and  data  type,  are  used. 

The  most  common  of  these  are  the  radar  sigmas.  The  residuals  and  partial 

derivatives  of  a  given  type  of  observation  from  a  specific  station  are  divided 

by  the  corresponding  <r  ..  In  this  way,  it  is  possible  to  insure  that  the  more 

SJ 

accurate  data  have  a  large  part  in  determining  the  optimum  set  of  parameters. 

If  it  is  known  (or  suspected)  that  a  station  has  a  constant  bias  in  reporting 

either  data  or  time,  the  P  .  are  used.  These  are  applied  to  the  appropriate 

SJ 

components  of  Omc  (before  the  division  by  Biases  may  be  included 

with  the  parameters  to  be  solved  for,  but  only  biases  on  basic  data  types. 

In  the  succeeding  paragraphs,  the  sigmas  and  biases  in  A  and  Omc  are  in¬ 
cluded  implicitly. 


3.  7.  3 


The  Unconstrained  Normal 


In  its  simplest  form,  differential  correction  involves  the  solution  of  the 

T  T  T 

linearized  problem  (A  A)  AP  =  A  O  .  A  A  is  the  normal  matrix.  If  a. 

r  T  me  l 

is  a.  row  of  A  and  if  A  A  =  0  initially,  the  normal  is  formed  in  TRACE  by 

T  T  T 

accumulating  A  A  =  A  A  +  a.  a^,  i  =  1,2,  ...»  n,  at  each  observation 
time . 
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3.  7.4 


The  Constrained  Normal 


It  is  often  desirable  to  impose  linear  constraints  of  the  form  AP  =  B{AP!)  +  C 
on  the  solution.  P'  is  Some  subset  of  P,  and  C  is  a  vector  of  constants.  For 
instance,  suppose  the  parameters  to  be  solved  for  are 


where  S ^ , 


are  two 


latitude 
longitude 
latitude 
longitude 
range  bias 

radar  stations. 


Then  the  requirement  that  the  positions  of  the  station  relative  to  each  other 
remain  fixed  is  equivalent  to  the  matrix  equation 


‘Apl " 

’l  0  0 

Apl 

Ap2 

0  1  0 

Ap2 

Ap3 

= 

1  0  0 

AP5 

-1- 

0 

Ap4 

0  1  0 

_  Ap5  _ 

0  0  1 

. 

AP  B  AP'  C  (65) 

T 

The  effective  problem  now  is  to  solve  the  reduced  system  (AB)  (AB)(AP')  = 

T  T 

(AB)  O  .  If  a.  is  a  row  of  A,  and  if  A  A  =  0  initially,  the  restricted  normal 

m  G  1  rp  rp  rp, 

is  formed  by  accumulating  A  A  =  A  A  +  (a.B)  (a.  B). 

3.7.5  Bounds 

Given  a  set  of  bounds  g^,  the  corrections  Ap.  to  the  components  of  P  are: 

(a)  less  (in  absolute  value)  than  g^  if  g^  >  0 

(b)  zero  if  g^  =  0 
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(c)  unrestricted  if  g^  <  0 
for  i  =  1,  2,  .  .  .  ,  m. 

If  constraints  are  to  be  applied,  the  bounds  are  adjusted  to  the  new  problem, 
m  (Sign  g  ) 

Leth  =;r  - - J—  b  j  =  1, 2 . k 

J  i=l  (g  /  1J 

Then  there  will  be  k  new  bounds,  g! ,  where 

g  •  _L—  if  h  >  0 

Vh  j 


g. '  =  0  if  h.  -  0 
J 


g.'  <  0  if  h.  <  0 
SJ  J 


Notice  that  g/  =  g.  for  a  variable  not  appearing  in  any  restraint 
equation.  Also,  specifying  bounds  that  are  equal  in  magnitude  but  op¬ 
posite  in  sign  for  two  parameters  to  be  corrected  by  equal  increments 
will  result  in  a  zero  correction  to  both. 


3.7.6  Solution  of  the  Normal  Equations 

For  the  purpose  of  this  section,  assume  that  there  are  m  para¬ 
meters  p^.p^,  ....  pm  to  be  corrected.  (The  succeeding  paragraphs 
require  no  change  other  than  a  substitution  of  k  for  m  in  the  constrained 
case.  ) 


T  T 

At  this  point  then,  we  have  the  mXm  matrix  A  A,  the  vector  A  Omc> 
and  a  set  of  bounds,  g_  The  problem  then  is  to  minimize  1 1 AAP  -  °mc  1 1  ^ 
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under  the  side  condition  that  h 


2 


for  which  g.  >0. 

°i 


the  sum  being  taken  over  all  i 


We  can,  without  loss  of  genei*ality,  assume  that  g.  4  0.  This  is  so  because 

th  ^ 

=  0  implies  that  Ap^  =  0,  and  the  i  row  and  column  of  the  normal  matrix 
may  be  ignored.  This  simply  reduces  the  dimension  of  the  problem. 


Now,  let  G  be  the  diagonal  matrix  so  that  G..  =  g^-  if  g.  >  0;  g. .  =  0  if  g.  <  0. 

We  wish  to  find  a  value  of  z  so  that  the  solution  AP'(z)  of  the  linear  system 
T  Z  T 

(A  A  +  zG  )  AP  =  A  O  satisfies  the  given  side  condition.  This  involves 

me  & 

two  procedures:  the  choice  of  the  best  value  for  z,  and  the  actual  solution 
of  the  system. 


3.7.6.  1  Determination  of  z 

As  a  start,  find  AP'(0),  (the  solution  to  (A^A)  AP  =  A^Omc).  If 


1  +  e  the  problem  is  solved.  (€.  and  €  here  are  suitably  small 


positive  constants.  )  If  not,  define  y(z)  =  Z 


/ Ap  (z)\ 

hr  J  -  u 


Now  y(0)  >  e  ^ 


Compute  y(h),  y(10h),  y(100h),  .  .  .  ;  h  some  present  constant,  until  either 


(1)  a  value  of  z(  =  kh)  is  found  so  that  <  y(z)  <  €]_,  in  which 
case  AP(z)  is  the  solution,  or 

(2)  two  values  of  z  are  found  so  that  y(z]_)  >  «]_  and  y(z^)  < 
The  required  value  of  z  is  now  bracketed. 

If  (2)  is  the  case,  the  next  step  is  to  choose  a  value  z^  between  z^  and  z^. 
z  ^  =  0.  8z  ^  -f  0.  2z^  whe  re  the  0.8,  0.  2,  are  fairly  arbitrary  and  z ^  and  z^ 
may  have  been  interchanged,  so  that  z^  is  closest  to  the  value  of  z  giving 
the  smallest  y(z). 


If  <  y(z^)  <  AP'(z^)  is  the  solution.  Otherwise,  inverse  quadratic 

interpolation  is  used  to  obtain  a  new  guess  z^.  Again,  if  <  y(z^)<€^, 

AP'(z.)  is  the  solution.  If  not,  the  uwo  values  of  z,  from  the  set  (z,,  z-,, 
4  1  c. 
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z y  z^),  which  bracket  the  solution  most  tightly,  are  chosen  and  the  process 
is  repeated  from  (2). 

If  more  than  twenty  solutions  of  the  linear  system  are  required,  the  process 
is  abandoned. 

3.  7.  6.  2  Solution  of  the  Linear  System 

This  section  describes  the  procedure  used  in  solving  the  linear  system 
T  2  T 

(A  A  +  zG  )  AP  =  A  O  .  For  a  discussion  of  the  theory  involved,  see 

me  7 

Section  2. 

T  2  T 

Let  C  =  A  A  +  zG  .  It  is  desired  to  find  a  matrix  S  so  that  SCS  =  D;  S 

lower  triangular  with  (-1)  on  the  diagonal;  D  =  diag  (d^,  ..  .,  d^). 

Il  C  is  a  one-by-one  matrix,  S  =  -1,  D  =  C.  Now  augment  C  by  another 
row  and  column: 


Since  S'  must  be  lower  triangular,  with  (-1)  on  the  diagonal,  and  of  the  same 
order  as  C',  it  must  be  of  the  form 


and  the  requirement  S'C'S’  =  D'  is  equivalent  to  solving  for  a  vector  W  and 
a  scalar  b  such  that 


3 -‘36 


It  is  easily  verified  that 


W=STD'1Sd  and 


b  =  a  -  W  d 


satisfy  the  requirements. 


The  computations  follow  the  above  outline,  starting  with  the  two-by-two 
matrix  C„,  i  =  1,  2,  ;  j  =  1,  2  and  continuing  until  the  decomposition 


T  2  T 
A-^A+zG  A  O 


WT  -1  /  \  O  1  A 
/  \  me 


T 

O  O 

me  me 


T 

Sl  W 


0  -1 


D  0 


0  a 


has  been  found. 

Carrying  out  the  multiplication  indicated  on  the  left  side  of  the  equation  shows 
that  the  m -dimensional  vector  \V  is  the  solution  to  the  linear  system. 

3.  7.  6.  3  Residuals  Prediction 

Hop  J  =  IIaap-o  II  is  computed  from  the  augmented  normal  matrix: 
me  me  ° 


AAP  -  O  1,2  =  f  (ATA  AP)  •  AP]  -  2  [<AT0  J  •  AP]  +  O  T  O  (68) 
mr  '  '  m  c  me  me 


3.  7.  6.  4  The  Inverse  Normal 


T  - 1  T  - 1 
(A  A)  1  =  S  D  LS 


with  S  and  D  as  defined  in  paragraph  3.  7.  6.  2. 


3.  7.  7 


Convergence  of  the  Differential  Correction  Process 


The  l!O  nc  [I  is  a  measure  of  how  well  the  orbit,  computed  on  the  basis  of  a 
given  set  of  parameters  P,  fits  the  observed  data.  Ho^^H,  computed  as 
in  paragraph  3.  7.  6.  3,  is  an  approximation  to  II  ^mr  II »  which  would  be  ob¬ 
tained  by  replacing  P  by  P  +  AP.  This  approximation  would  be  exact  if  the 
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least  squares  problem  were  linear;  that  is,  if  P  were  in  a  sufficiently  small 
neighborhood  of  the  minimum  point. 


Convergence,  then,  is  defined  as  being  that  point  at  which  further  corrections 
corrections  to  P  would  produce  no  significant  decrease  in  ||  0^c  [! ;  i.e.  ,  no 
over-all  improvement  of  the  residuals.  The  criteria  used  are 


O  ||  -  ||  op  || 


me 


me 


O 


<  c  .  or  O 
—  1  me 


•  n  where  n  is  the  number  of 


me 


observations  and  e  ^  and  are  input  quantities 


If  i]  O  c  ||  is  decreasing  with  each  iteration,  the  process  is  converging 
and  the  bounds  are  expanded  at  each  step  (by  a  multiplicative  factor  (3^)  to 
permit  faster  convergence.  On  the  other  hand,  if  ll^mcll  is  increasing  from 
one  iteration  to  the  next,  the  process  is  diverging  and  the  last  corrections  are 
presumed  to  have  altered  P  too  drastically.  In  this  situation  the  previous 
values  of  P  and  the  corresponding  normal  matrix  are  retrieved  and  resolved 
with  tighter  bounds.  The  new  bounds  g!  are  such  that  the  weighted  length 
"G1  •  AP1  of  the  solution  is  reduced  to  6^  times  its  previous  value. 

3.7.8  The  Correlation  Matrix 

If  the  mathematical  model  is  exact,  if  the  observations  are  linear 
functions  of  the  parameters,  if  the  observation  errors  have  mean  zero 
and  are  independent,  and  if  the  input  values  of  are  correct,  then  the 
inverse  of  the  normal  matrix  is  the  variance-covariance  matrix  of  the 
parameters,  due  to  the  random  errors  in  the  observations.  (For  a  proof 
of  this,  see  Section  2.  )  If  the  elements  of  this  matrix  are  given  as  c^, 
the  corresponding  correlation  matrix  has  elements 


c  ? . 
ij 


c. . 

U 


c . . 
JJ 


j. 


1,  2, 


m. 


(70) 
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If  all  of  the  above  assumptions  are  true  except  that  all  of  the  tr^ 
values  are  in  error  by  a  constant  multiplicative  factor,  then  the  values 
in  the  variance-covariance  matrix  will  also  all  be  in  error  by  the  same 
factor. 
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3.8 


ERROR  ANALYSIS 


It  is  often  desirable  to  analyze  the  basic  statistics  involved  in  a  particular 
orbit  determination  problem.  This  essentially  entails  a  determination  of 
the  effects  that  specified  sources  of  error  have  on  the  precision  of  the  least 
squares  parameters.  These  sources  of  error,  for  example,  may  be  in¬ 
accuracies  in  station  locations,  random  errors  in  observations,  or  errors 
in  differential  equation  parameters. 


The  error  analysis  procedure  does  not  require  an  actual  determination  of 
the  orbit,  nor  does  it  require  actual  observations;  however,  the  data  types 
and  data  intervals  must  be  specified  for  each  station.  A  basic  error  analysis 
may  then  be  obtained  by  simple  matrix  manipulation. 

3.  8.  1  Notation  and  Nomenclature 


n  number  of  "observed"  quantities 

m  number  of  least  squares  parameters 

k  number  of  parameters  (other  than  least  squares  parameters), 

which  are  considered  in  error 

cr  .  radar  sigma  to  be  applied  to  data  type  j  from  station  s 
s) 

diagonal  matrix  of  order  n  with  <r”  }  for  each  "observation" 

&  sj 

as  elements 

P  vector  of  least  squares  parameters,  or  "P-parameters" 

Q  vector  of  parameters  ("Q-parameters")  considered  in  error 

(other  than  least  squares  parameters) 

C(Q)  variance- covariance  matrix  of  the  specified  Q-parameters 
Ap  n  X  m  matrix  of  weighted  partial  derivatives  of  the  observa¬ 
tions  with  respect  to  the  least  squares  parameters; 

1/2  OOc 
A  =  W1^ 
p  oP 

A  n  X  k  matrix  of  weighted  partial  derivatives  of  observations 
with  respect  to  Q-parameters 


A 
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3.  8.  2 


The  Normal  Matrix 


T 

The  normal  matrix  is  formed  in  TRACE  by  accumulating  A  A  = 

^  T  th 

a.  a.  where  a.  is  the  i  row  of  the  matrix  A  =  (A  ,  A  ).  In  terms 

i=i  1  1  1  p  q 


of  P  and  O  parameters  the  A  A  matrix  is: 


(71) 


This  accumulation  is  made  in  double  precision  in  error  analysis  applications 
only. 


3.  8.  3  (A  TA  )  1 

P  P 


The  matrix  is  the  variance-covariance  matrix  of  the  least  squares 
(P)  parameters  due  to  random  errors  in  the  observations.  (The  random 
errors  in  the  observations  are  specified  by  the  o,gj  .  )  In  the  case  of  the 
orbital  parameters,  the  uncertainties  given  by  this  covariance  matrix 
C(P'}  =  (A*p^A^)  ^  apply  only  at  the  initial  time  tQ. 
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3.  8.  4 


T  - 1 

Effect  of  Q-parameter  Errors  on  (A  A  ) 

_ _  _  P  P 

In  paragraph  3.  8.  3,  only  the  effect  of  observation  errors  on  the  set 
of  least  squares  parameters  was  considered.  However,  station  location 
errors  and  differential  equation  parameter  errors  will  also  contribute  to 
the  uncertainty  of  the  P-par amete rs.  By  including  these  ''Q"  effects  a 
new  covariance  matrix  is  given  by: 


C(P")  =  (A  TA  )_1 

P  P 


T  - 1 
+  (A  1 A  )  1 
P  P 


A  TA  C(Q.)A  TA 

P  q  q  p 


T 

(A  1 A  ) 
P  P 


-1 


=  C  (P ' )  + 


(72) 


where 


dP" 


T  - 1 
=  -(A  'A  )  1 
P  P 


Again,  this  variance-covariance  matrix 


applies  at  the  initial  time  tQ.  (Note:  The  C(Q)  matrix  contains,  essentially, 
the  uncertainties  in  the  O-parameters  and  must  be  input.  ) 


3.8.5  Transformation 

Both  C(P')  and  C(P")  are  referred  to  the  initial  time.  Since  un¬ 
certainties  in  initial  coordinates  do  not  satisfactorily  describe  trajectory 
uncertainties,  it  may  be  desired  to  translate  these  covariance  matrices 
into  other  coordinate  systems  and  to  times  other  than  tQ.  This  is  some¬ 
times  called  the  "updating"  process.  The  general  transformation  equation 
is 


C(Xt) 


ax  aP" 
ap  ao 

o  o 


C  (Q) 


ax 

^"o 


ax  aP"  \T 
ap  ao  ) 

o  o  / 
(73) 


Transformations  to  other  coordinate  systems  are  accomplished  as  follows. 
(In  all  cases  the  result  is  a  sum  of  two  terms:  The  first  gives  the  effect 
of  observational  errors  only;  O-parameter  effects  arise  from  the  second. 
Both  the  first  term  and  the  sum  can  be  printed.  ) 
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C(Rt) 

C(Et) 

C(Tt) 


(orbit  plane) 


(spherical) 


(element) 


(period,  apogee,  perigee) 


(74) 

(75) 

(76) 

(77) 


3.8.6  Transformation  Partial  Derivatives 

The  following  formulae  are  used  for  the  transformation  partial  derivative 
matrices  in  paragraph  3.  8.  5. 

Orbit  plane  coordinates: 

%  =  -  (where  X  =  (  y  )  for  this  and  the  following  equation  only) 

I  X  |  \  z  / 


gxX 

U*x| 
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b 


X 
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Let  ’ii  be  the  matrix  ,  T|,  Q  ),  then 


Spherical  coordinates: 
a  (right  ascension) 


do-  _  -y 

dx  '  2  2 

x  -  y 

bo  x 


bo  x 


6  (declination) 


i*6  -  xz 


bt  _  -  yz 


*L'rt 


j3  (flight  path  angle) 


SB  1 
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SB 

sx 


-  X 


A  (azimuth) 


SA 
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y  (rz-zr)  -  (xy-yx)  (xz-  zx  +  ) 

7  2  .2..  2  27 

(v  -r  )  (x  +y  ) 


SA 

Sy 


-x  (rz-zr)  +  (xy-yx)(yz-zy  +  X^£)  i 

2  2  2  2. 

(v  -r  )  (x  +y  ) 


SA  _ 

r  (xy-yx) 

Sz 

2.  2  .2. 
r  (v  -r  ) 

SA  _ 

(yz-zy) 

Sx 

2  .2. 

r^v  -r  ) 

SA  _ 

(xz - zx) 

,  2  .2. 
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:r  (  SA  ' 

Sz 

r  y  Sz 

r  (radius) 


v  (velocity) 


a  (semi-major  axis) 


da 

dr 
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2 

where  X.  = - 
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where  p  =  a  (1  -e  ) 
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=  sin  (a-Q) 


Si 

SA  “ 


-  cos  5  cos  (o'-O) 


Si  _  _Si_  _  _Si_  _  Si_ 
Scy  ~  S3  Sr  -  ,Sv 


0  (longitude  of  ascending  node) 


SO  _  -cos  (o'-O) 
SS  ~  tan  i 


SO  _ 
SA 


-sin  6 


.  2  . 
sin  i 


so 

So- 


1 


so  ^  so  _  so 

S3  Sr  Sv 


0)  (argument  of  perigee) 


Suj  _  cos  (cv-Q) 
S6  cos  A 


-N  X 

Su  _  OJ 

SB  "  ae 


A  t 


where  x 

U) 


p-r 


e 


Sm  _  cos  6  sin  (o'-O ) 
SA  “  sin  i 

Su;  y« 


3-67 


bv  Z2u_ 
dv  e  rv 


t  (time  of  perigee  passage) 


St  -\/P  Xu) 
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e  yja 


dT 
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dE 
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where  M  is  mean 
anomaly 
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(aP°gee) 


dr  r 

a  a  , 

55-  =ir  =  1+e 


Sr. 

de~ 


=  a 


rp  (perigee) 


dr  r 

T5-E  =  -E  =  I-e 

da  a 


dr 


de 


B  =  -a 


dT 


(All  other  terms  of  are  zero.  ) 
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SECTION  4 


PROGRAM  STRUCTURE 


4.  1  GENERAL 

TRACE  is  written  in  the  FORTRAN  language,  to  be  used  with  the  IBM 
709/7090  FORTRAN  Monitor  System.  The  ba  sic  construction  of  the  program 
is  a  series  of  independent  links,  which  are  connected  by  the  CHAIN  feature 
of  FORTRAN.  Within  each  link  there  is  a  series  of  large  blocks,  or  major 
subroutines,  each  of  which  makes  use  of  many  smaller  routines. 

This  design  was  chosen  for  TRACE  for  several  reasons.  First,  the  flow  of 
computation  is  easy  to  follow  and  understand,  both  in  general  and  in  detail, 
by  relative  newcomers  as  well  as  by  the  authors.  This  is  an  important 
consideration  in  facilitating  modifications  to  an  intricate  but  continually  ex¬ 
panding  program  such  as  TRACE.  Because  of  the  subroutine  structure  of 
TRACE,  most  of  the  presently  projected  additions  can  be  made  by  changing 
one  or  two  routines  with  no  possibility  of  interfering  with  contiguous 
segments. 

TRACE  is  restricted  to  use  with  the  709/7090  by  only  one  characteristic — 
the  occasional  utilization  of  FAP.  In  general,  the  FAP  subroutines  are  short 
and  few.  In  only  one  link  (TRAIN,  the  tracking  data  input  link)  do  they  play 
a  major  part,  and  even  in  this  case  they  can  be  simply  replaced.  It  was  felt 
desirable  to  replace  the  FORTRAN  input  statements  by  buffering  routines 
designed  specially  to  handle  large  amounts  of  fixed  format  radar  observa¬ 
tion  data  cards. 

TRACE  is  therefore  an  extremely  flexible  program  partitioned  into  five 
major  links.  A  description  of  these  links  follows  and  general  flow  charts 
(Figures  10  through  15)  appear  at  the  end  of  this  section. 
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4.  2 


THE  LINKS 


4.  2.  1  CHAIN 

This  is  the  only  program  that  must  be  executed  regardless  of  the  mode  in 
which  TRACE  is  to  be  used,  CHAIN  is  divided  into  three  sections.  The 
first  section  reads  basic  data,  but  not  station  locations,  or  tracking  data, 
or  the  specifications  for  data  generation  or  error  analysis.  It  also  prints 
a  header,  sets  several  options  to  their  nominal  values,  and  computes  the 
Julian  Date  and  the  orientation  of  the  Earth.  The  second  section  initializes 
the  trajectory  parameters,  and  reads  in  extra  input  for  GAIN  or  FEIGN, 
if  the  corresponding  mode  is  being  executed. 

The  third  section  is  executed  only  during  the  orbit  determination  mode  of 
MAIN  {see  paragraph  4.  2.  3)  .  It  contains  the  differential  correction  pro¬ 
cedure  and  the  corresponding  output.  Control  is  alternated  between  MAIN 
and  this  section  of  CHAIN  during  the  orbit  determination  mode. 

4. 2.  2  TRAIN 

The  tracking  input  link,  TRAIN,  reads  radar  station  location  and  observa¬ 
tion  data.  This  data  may  be  on  the  BCD  input  tape  produced  by  the  IBM  1401, 
or  on  a  binary  tape  previously  written  by  TRAIN.  For  real-time  tracking 
exercises,  the  card  reader  will  be  used.  The  observations  are  sorted 
chronologically,  and  a  compacted  list  is  produced,  which  eliminates  storage 
corresponding  to  blanks  in  the  information  reported.  In  this  way,  approxi¬ 
mately  two  thousand  observations  can  be  handled  in  core  {on  a  32K  machine) 
without  resorting  to  intermediate  tapes. 

On  option,  TRAIN  produces  a  binary  tape  containing  the  sorted  and 
compacted  observation  data,  to  be  used  either  on  successive  runs  or  the 
next  case  of  the  same  run.  This  data  has  been  transformed  to  the  standard 
set  {  R,  A,  E,  R,  P,  Q,  P,  O  at  present),  and  all  units  have  been  converted 
to  an  earth-radii-minutes-radians  system.  Observation  times  are  reduced 
to  minutes  from  midnight  of  epoch. 
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In  addition,  TRAIN  prints  the  tracking  input,  and  decodes  and  prints  input 
concerning  parameters  to  be  solved  for  by  differential,  correction. 

4.  Z.  3  MAIN 

MAIN  has  two  modes  of  operation:  trajectory  only,  and  orbit  determination. 
The  first  is  a  straightforward  compulation  of  the  trajectory  determined  by 
the  initial  conditions,  using  whichever  of  the  available  methods  has  been 
designated.  Output  is  a  time  history  of  position  and  velocity  of  the  vehicle, 
and  various  other  related  quantities,  at  any  specified  combination  of  print 
intervals. 

The  orbit  determination  mode  utilizes  both  the  trajectory  block  and  a  radar 
block.  These  two  segments  combine  to  produce  residuals  and  partial 
derivatives  of  the  observations  with  respect  to  the  parameters  being  studied. 
The  original  nonlinear  problem  is  solved  iteratively  by  differential  correc¬ 
tion.  Residuals  are  formed,  corrections  computed  and  applied  and  the  entire 
process  is  repeated  either  a  specified  number  of  times  or  until  convergence. 
Convergence  is  defined  to  be  the  point  at  which  no  further  improvement  can 
be  predicted  in  the  residuals  (the  differences  between  measured  and  computed 
observations).  Divergence,  defined  as  an  over-all  increase  in  the  residuals, 
is  also  possible.  In  this  case,  the  last  best  solution  is  retrieved  and  the 
corresponding  system  re -solved  with  more  stringent  bounds  on  the 
corrections . 

Output  from  MAIN  in  this  mode  consists  of  the  residual  at  each  observation 
time,  and  the  pertinent  results  from  the  differential  correction  routines 
including  initial  conditions,  corrections,  and  the  correlation  matrix.  The 
trajectory  output  and  a  time  history  of  the  normal  matrix  and  its  inverse  are 
also  available  if  desired. 

4.2.4  GAIN 

Steering  data  for  the  radar  stations  must  be  supplied  in  the  form  of  separate 
listings  for  each  station.  Although  the  computations  involved  are  basically 
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those  of  the  computed  observations  produced  chronologically  in  MAIN,  the 
required  output  form,  station-by-station  listings,  imposes  a  considerable 
sorting  problem. 

Utilizing  the  trajectory  block  and  a  simplified  radar  block,  GAIN  proceeds 
so  as  to  completely  fill  core  with  chronological  ephemeris  data,  which  is  then 
sorted  into  listings  that  are  station-by-station  for  the  period  of  the  data  in 
core.  (Further  sorting  could  be  added  upon  completion  if  a  station-by- 
station  listing  for  a  longer  time  period  is  desired.  ) 

GAIN  requires  as  input  a  list  of  station  locations  and  a  definition  of  the  data 
configuration  desired  from  each.  Any  function,  not  necessarily  an  observa¬ 
tion,  could  be  coded  and  selected  for  output.  Normally,  however,  output 
will  be  some  type  of  radar  observation.  Data  is  produced  only  for  those 
periods  during  which  the  vehicle  is  visible  to  the  station;  optionally,  rise 
and  set  times  only  are  calculated. 

4.2.5  FEIGN 

The  simulation  link,  FEIGN,  is  designed  to  permit  studies  of  the  large 
matrices  involved  in  tracking  system  design  without  requiring  their  genera¬ 
tion  by  one  program  and  manipulation  by  another.  To  save  storage  space, 
the  simulated  data  is  not  computed  explicitly,  but  is  instead  inferred  from 
a  list  of  data  types  and  frequencies  for  given  stations.  FEIGN  employs  the 
trajectory  and  radar  blocks  exactly  as  they  are  used  in  the  differential 
correction  path  of  MAIN.  However,  since  there  are  no  actual  observations, 
this  link  checks  for  visibility  at  each  point  and  computes  derivatives,  etc.  , 
only  when  physically  applicable. 

The  main  matrix  calculation  being  made  in  FEIGN  at  present  is  the  calcula¬ 
tion  and  inversion  of  the  normal  matrix  associated  with  the  parameters  being 
studied  (i.  e.  ,  P-parameters ) .  The  inverse  of  the  normal  matrix  is  the 
covariance  matrix  of  the  P-parameters.  This  covariance  matrix  may  be 
updated  in  time  (by  use  of  the  variational  equation  partials)  and  transformed 
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to  other  systems.  These  other  systems  are:  Cartesian;  Orbit  Plane; 
Spherical;  Element;  and  Period,  Apogee,  Perigee.  The  effect  of  specified 
parameter  errors  (Q-parameters)  may  be  included  in  any  of  the  above 
matrices.  Any  combination  of  matrices,  with  or  without  the  Q-parameter 
effects,  may  be  output. 
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Figure  10.  General  Flov/  Chart  for  TRACE 


b 
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CHAIN 


Figure  11.  Flow  Chart  of  CHAIN  Link 

t 
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TRAIN 


Figure  12.  Flow  Chart  of  TRAIN  Link 
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MAIN 


Figure  13.  Flow  Chart  of  MAIN  Link 
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GAIN 
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FEIGN 


Figure  15.  Flow  Chart  of  FEIGN  Link 


4-11 


FEIGN 


Figure  15.  Flow  Chart  of  FEIGN  Link 
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SECTION  5 

USAGE 

5.  1  INTRODUCTION 

This  section  describes  the  input,  deck  setup,  and  output  connected  with  the 
use  of  TRACE  (in  Sections  5.2,  5.2.6  and  5.3,  respectively).  It  is  designed  to 
answer  most  questions  pertaining  to  the  use  of  the  program.  Each  function- 
trajectory,  tracking,  data  generation,  and  error  analysis  -  is  treated 
separately.  A  detailed  explanation  of  each  is  contained  in  Section  1. 
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5.2 


INPUT 


This  section  is  divided  into  six  parts.  The  first  is  concerned  with  the 
Basic  Data  input  common  to  all  functions.  The  next  four  parts  cover 
the  functions  of  trajectory,  tracking,  data  generation,  and  error  analysis. 
Each  part,  together  with  the  Basic  Data  input,  is  independent  and  completely 
describes  the  input  for  the  function  involved.  Finally,  the  arrangement  of 
the  program  and  input  deck  is  described. 

TRACE  utilizes  the  following  four  types  of  load  sheets  for  input: 

•  FINP  for  Basic  Data 

•  Station  Location  Data 

•  Radar  Observation  Data 

•  Error  Analysis  and  Data  Generation  Specification 

Sample  load  sheets  and  summary  descriptions  are  included  at  the  end  of 
each  part. 


The  following  three  points  of  information  will  make  the  FINP  load  sheets 
easie r  to  use. 

a.  Although  the  load  sheet  imposes  an  order  on  the  input,  the 
actual  order  of  the  cards  is  (almost)  immaterial,  the  only 
restriction  being  that  all  values  with  nonsymbolic  locations 
are  located  relative  to  the  last  previous  symbol.  In  the 
case  of  two  appearances  of  the  same  location  symbol,  the  last 
value  read  is  the  effective  input. 

b.  A  prefix  (Columns  1,  19,  37,  and  55)  determines  the  mode 
of  input.  A  blank  indicates  that  the  following  value  is  to  be 
read  as  a  floating  point  number;  an  I,  as  a  fixed  point 
decimal;  D,  BCD  (Hollerith);  and  B,  as  an  octal.  The  END 
cards  are  used  because  the  prefix  E  terminates  the  FINP 
read. 

c.  Any  card  for  which  no  value  appears  may  be  omitted.  Blank 
fields  are  ignored  except  for  D  prefix  (BCD). 
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5.  2.  1 


Basic  Data 


By  definition,  Basic  Data  is  that  which  is  common  to  all  functions.  The  re¬ 
quired  Basic  Data  includes  the  list  of  functions  to  be  performed,  the  specifi¬ 
cation  of  the  trajectory  (date,  time,  and  initial  conditions},  the  force  model 
to  be  assumed,  and  the  constants  and  parameters  to  be  used  in  ihe  trajectory 
integration.  (The  force  model,  constants,  and  parameters  are  "required" 
data,  but  standard  values  are  provided.  The  replacement  of  these  quantities 
is  thus  'optional."  See  the  appendix  for  a  list  of  the  standard  values.) 

There  are  also  options,  which  are  common  to  all  functions,  and  thus  are 
contained  in  the  Basic  Data  input.  These  include  identification  information, 
specification  of  the  ballistic  (drag)  coefficient  and  atmosphere  model,  se¬ 
lection  of  other-body  perturbations,  and  output  print  time  specification. 


A  line-by-line  explanation  of  the  load  sheet  follows.  An  asterisk  refers  to 
a  special  feature,  which  is  explained  below  the  example. 

5.  2.  I.  1  Required  Input 

Line  1  -  Functions  To  Be  Performed 


1  UD  \  1TIN  1 


2  4 


This  line  contains  the  ordered  list  of  all  functions  TRACE  is  expected  to 
perform  on  a  given  run,  using  the  code: 


tracking  data  input  =  1 

tracking  computations  =  2 

trajectory  only  =  3 

data  generation  =  4 

error  analysis  =  5 


Up  to  12  functions  may  be  selected.  When  this  list  is  exhausted,  TRACE 
will  reset  certain  standard  options,  prepare  to  run  another  sequence  of  func¬ 
tions,  and  read  Basic  Data,  or  stop  if  there  is  none.  The  example  given  would 
be  for  a  tracking  case  (functions  1  and  2)  followed  by  a  data  generation  case. 
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Lines  2  through  8  -  Epoch 


*  2 

1 

YEAR 

1963 

3 

I 

MNTH 

6 

4 

I 

DAY 

15 

5 

TZNE 

0 

6 

HR 

12 

t 

MIN 

45 

8 

SEC 

15.  5 

The  year,  month,  and  day  should  be  input.  The  X  axis  is  then  directed  to  the 
vernal  equinox  (see  paragraph  3.  1.  1).  (*  Line  2.  If  the  year  is  input  nega¬ 

tive,  the  X  axis  will  be  directed  to  the  longitude  of  Greenwich.)  The  hour, 
minute,  and  second  refer  to  midnight,  zone  time.  GMT  is  time  zone  zero. 


Lines  9  through  15  -  Initial  Conditions 


9 

1 

ICTYP 

2 

10 

IC 

126.  1 

1 1 

2 

31,  23 

12 

3 

89- 

1  3 

4 

14. 

*  14 

5 

22600114 

*  15 

6 _ 

25117.  3 

Line  9  indicates  the  type  of  initial  conditions  entered  in  lines  10  thru  15. 
For  ICTYP  =  :  IC's  are: 


1  Earth-centered  inertial  cartesian  coordinates  (x,  y,  z,  x,  y,  z 
in  feet  and  ft/sec  -  see  paragraph  3.  1.  1.  1). 

2  Spherical  coordinates  (a,  6,  p,  A,  r,  v  in  degrees,  feet  and 
ft/ sec  -  see  paragraph  3.  1.  1.2). 

(:‘:Line  14.  If  r  is  negative,  it  is  interpreted  as  height 
above  the  earth's  surface.  ;!'Line  15.  If  v  is  negative, 
circular  velocity  is  computed  and  used.  ) 

3  Orbital  elements  (a,  e,  i,  ft, «,  t  in  feet,  degrees,  and 
minutes  -  see  paragraph  3.  1.  1.  3). 

4  The  same  as  2  above  with  longitude  replacing  right 
ascension. 
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5  No  IC's  are  input.  The  last  trajectory  point  of  the  last 
previous  run  is  used. 

6  No  IC's  are  input.  The  corrected  initial  conditions  from 
the  last  previous  tracking  run  are  used. 

8  Same  as  1,2,  and  3  but  in  units  of  earth  radii,  minutes, 
and  radians  The  type  is  determined  from  the  last  previous 
tracking  run,  or  from  CPRAM  (see  line  37). 

9  Same  as  1  but  in  units  of  earth  radii  and  earth  radii/min. 

0  No  IC's  are  input.  .For  a  tracking  run  two  RAE  sets  are 
used  from  the  data  to  calculate  a  set  of  initial  conditions 
(see  paragraph  3.6.7). 

5.  2.  1 .  2  Optional  Input 

Lines  16  through  20  -  Drag 


1 6 

DRAG 

.  01 

17 

I 

2 

18 

3 

19 

4 

20 

5 

D  2 

Line  16  contains  the  drag  parameter — in  ft  /lb. 


Line  17  contains  the  atmosphere  model  specification. 

If  DRAG(2)  is:  0  use:  ARDC  59  Model 

1  Lockheed  Model 

2  Paetzold  Model  II 

3  Paetzold  Model  I 

4  L.  F.  E.  Model 

If  line  17  is  a  2,  3,  or  4,  lines  18  through  20  contain  three  quantities  used 
in  the  Paetzold  calculation. 


Line 

Line 

Line 


18  contains 

19  contains 

20  contains 


F  -  solar  radio  flux. 

Ap  -  planetary  magnetic  index. 
g(a)  -  plasma  intensity  coefficient. 
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Drag  Table  Option 


cnA  cnA 

Think  of  the  drag  parameter  as  two  terms:  ^  X  C^(M,  h).  — is  a 

constant  and  can  be  differentially  corrected  by  the  use  of  the  variational  eq¬ 
uation.  C^(M,  h)  is  a  function  of  Mach  No.  and/or  altitude;  it  is  then  input  as 
a  table . 


Normally  Cj^  is  set  equal  to  1,  and  with  — input  into  location  DRAG, 
TRACE  operates  as  before.  If  the  use  of  a  table  is  desired,  additional  inputs 
must  be  made  {marked  with  *  below).  The  following  is  the  use  of  the  C  block 
where  the  tables  are  stored. 


* 


C(50) 

C(51) 


^(52) 


C(53) 

C(54) 

C(55)  -  C(7Z) 

C(73) 

C(74) 

C( 75)  -  C(100) 


Cj^(=l,  or  interpolated  table  value) 

if  =  0,  do  not  use  tables 

if  negative,  use  Mach  table  only 

if  positive,  use  Mach  table  and  altitude  table 

altitude  above  which  altitude  table  is  used  and 

below  which  Mach  table  is  used.  (Needed  if 

C(5  1)  is  positive .  ) 

(not  used) 

(used  by  interpolation  routine) 

altitude  table.  (Monotonic  increasing  stored  as 
altr  C^,  alt.,,  CD  .  .  .  altn,  C^,  0,  0) 

(not  used) 

(used  by  interpolation  routine) 

Mach  No.  table  (stored  as  altitude  table  above). 


Line  21  -  Print  Code 


■ 

■ 

is  cu<  >  <u  u  d NT) 

21 

_ 

D 

PRCDE 

X 

■ 

■■ 

S3 

■Ql 

■ 

■ 

■ 

N / 

■5 

SIS 
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An  X  will  cause  the  corresponding  information  to  be  output.  (The  boxes 
marked  V  will  be  printed  only  at  specified  print  times  -  see  below.  }  See 
Section  5.3  for  output  examples. 

Lines  22  through  28  -  Print  Time  Vector 


22  1 

I 

PRTIM 

0 

i 

n 

23 

I 

2 

2 

24 

3 

0 

Ati 

25 

4 

5 

tl 

26 

5 

90 

*  *2 

27 

6 

15 

zZ 

28 

7 

1440 

Sequence  of  print  times  for  output  selected  in  PRCDE.  There  may  be 
n{n  <  9)  sets  {line  23);  for  the  ith  set  output  is  from  t.  ,  to  t.  at  intervals  of 
At..  All  times  are  in  minutes  from  midnight  if  PRTIM  =  1;  from  epoch  if 
PRTIM  =  0.  At.  =  C  means  do  not  print  in  this  interval.  Additional  cards 

i  r 

may  be  inserted  here  if  3  <  n  <  9. 


Line  33  -  Extra  Body  Perturbation 


33  [  I 

CTAPE 

7 

If  CTAPE  is  non-zero,  then  other-body  perturbations  will  be  computed 
from  coordinate  tape  on  unit  number  CTAPE,  using  body  selectors  in  INTEG 
and  relative  masses  in  C. 


Line  34  -  Eclipse  Indication 


3~4BT[YfAPE 


If  XT  APE  is  non-zero,  then  the  position  of  the  sun  will  be  determined  at 
each  print  time  from  the  coordinate  tape  on  unit  number  XT  APE.  An  indi¬ 
cation  is  then  printed  as  to  whether  the  sun  is  visible  from  the  satellite. 
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5.  2.  2 


Trajectory  Only 


5.  2.  2.  1  Required  Input 

There  is  no  required  input  in  addition  to  the  Basic  Data. 
5.  2.  2.  2  Optional  Input 


Lines  37  through  39  -  Variational  Equation  Partial  Derivatives 


An  X  in  a  box  causes  the  variational  equation  for  the  corresponding  parameter 
to  be  solved.  The  partial  derivatives  are  printed  out  if  the  Variational  Equa¬ 
tion  box  is  marked  in  PRCDE  (see  Line  21). 


CPRAM  -  Initial  condition  parameters.  The  first  box  specifies  the  type  of 
initial  conditions.  The  succeeding  boxes  indicate  the  particular  parameters 
desired.  The  boxes  are  ordered  as  follows: 


DPRAM  and  OPRAM  -  Differential  equation  parameters.  The  boxes  are 
ordered  as  follows: 


❖ 


DPRAM  Drag 

(L 

J2 

J3 

J4 

J5 

J  2 1  | 

J31 

J  . .  ! 
<*1 

J22 

J32 

J42 

OPRAM  J„ 

J43 

J44 

X21 

X31 

X41 

x22 

X32 

X42 

X33 

X43 

X44 

* 


This  list  applies  only  if  the  full  set  is  used.  See  the  appendix  for  the  ordering 
of  shorter  sets. 
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Line  42  -  Trajectory  Comparison  Option 


41 


42 


n 

IFLAG  1 

T] 

15 

10 

If  IFLAG{15)  is  +  it  contains  the  logical  tape  number  of  the 

tape,  which  is  to  contain  the  reference 
trajectory 

0  a  regular  trajectory  run  is  indicated 

it  contains  the  logical  tape  number  of  the  tape 
containing  the  reference  trajectory.  The  ref¬ 
erence  trajectory  is  read  in  and  differenced 
with  the  trajectory  of  the  present  case.  These 
differences  are  resolved  into  the  orbit  plane 
system  and  are  placed  on  the  binary  plot  tape. 
The  binary  plot  tape  number  is  IFLAG(15)  +  1. 

Line  44  -  Analytic  Trajectory  Option 


43 

1 

C 

P 

44 

U 

49 

LJ 

If  C(49)  is  non-zero,  the  trajectory  is  not  integrated.  Instead,  analytic 
approximation  formulae  are  used  to  determine  the  trajectory  (see  paragraph 
3.  6.  6).  C{49)  should  be  a  positive  integer;  then  the  formulae  are  recomputed 

every  n^  orbit  where  n  =  C(49). 


i 
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Table  1.  TRACE  -  Basic  Data  and  Trajectory  Input 


Table  2.  Basic  Data  and  Trajectory  Input  Description 


l:  Itinerary.  Specify  order  of  compu¬ 

tations. 

1 -Tracking  data  input,  2-Track, 

3- Trajectory  only, 

4- Data  generation, 

5- Error  analysis 

2-8:  Epoch 

9:  Type  1  2  3  4  of  initial  conditions 

10-15:  1C  x  o  a  ? 

y  6  e  6 

z  B  i  0  (ft-sec-deg) 
x  A  0  A 
y  r  j)  r 
z  v  t  v 

NOTE:  If  r  is  negative  it  is  inter¬ 
preted  as  height;  if  v  is  neg.  ,  cir¬ 
cular  velocity  is  computed  and  used. 
16:  CdA/W.  (ft^/lb) 

17:  Atmosphere  model  specification 

0  ARDC  59 

1  Lockheed 

2  Pactzold  II 

3  Paetzold  I 

4  L.  F.E. 


18: 

F 

-  solar  radio  flux 

19: 

A 

p 

-  planetary  magnetic  index 

20: 

g(a) 

-  plasma  intensity  coeff. 

§ 


21:  X  to  print 

22-28:  Print  at  t=t0(Atj)  tj  (At2).  .  .  <£tn)tn  in 
minutes  from  midnight  if  PRTIM-J  , 
from  epoch  if  0. 


t.  «  u 

2  rt  “  u 

.2,2 

2  «  z 


■  2  2  o  « 

<T  c  C  O'  G 

a)  5  m  u  (t  .13 

•  g  w  rt  T. 

11  S  coo  «i 

n!  ^  0  0-  li  .O 


i-  o.<3  >  o  o  a  n 


Initial  Condition  Parameter  Specification: 
TYPE  [5T]  to  specify: 

1  xyzxyzt 

2  o  6  0  A  r  v  t° 

3aei  Cl  u  t  t° 


38-39:  Differential  Eq.  Parameter  Specification: 

X  to  specify: 

cda  * 

W  U  J2  J  3  J4  J 5  J21  J31 


DPRAM 


J41  J22  J32  J42 


J 33  J 43  J44  L21  L31  L41 


OPRAM 


L22  L32  L42  L33  L43  L44 


42:  Trajectory  comparison  option  and  tape 

specification. 

44:  Analytic  trajectory  option. 

*  This  list  applies  if  <t>BJT  =  4.  If  <J)BJT  =  3 
replace  with  shorter  list  J 2J .  J-jj,  ^32’  ^33 

and  L-,^,  ^31’  ^22’  L32'  ^33'  ^EJT  =  2 
replace  with  J^,  J 22  and  L^  j ,  L^  2 . 


23  -  n  (<9),  24  -  t  ,  25  -  At,  . 

—  o  i 

26  -  tj,  27  -  At2,  28  -  t2>  etc. 

33:  Extra  body  perturbation  option 

and  coordinate  tape  number. 

34:  Eclipse  indication  option 

and  coordinate  tape  number. 
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5.2.  3 


T  racking 
5.  2.  3.  1  Requi red  Input 

Lines  33  through  37  -  Parameter  Specification  Boxes 


33 

4 

ia 

n 

a 

B 

B 

3 

■ 

— 

k  34 

13 

DPRAM 

■ 

■ 

■ 

s 

3 

3  5 

H 

CbPRAM. 

■ 

3 

B 

B 

■ 

1 

■ 

■ 

■ 

36 

E3 

12 

3 

B 

_ 

37 

_D 

m 

a 

Lj 

I 

_ 

- 

: 

!l 

J 

An  X  in  a  box  causes  the  corresponding  parameter  to  be  used  in  the  differen¬ 
tial  correction  solution, 

CPRAM  -  Initial  condition  parameters.  The  first  box  specifies  the  type  of 
initial  condition.  The  succeeding  boxes  indicate  the  particular  parameter 
desired.  The  boxes  are  ordered  as  follows: 


Type 


B 

II 

B 

2 

6 

8 

A 

3 

1  a 

e 

i 

B 

DPRAM  and  OPRAM  -  Differential  equation  parameters.  The  boxes  are 
ordered  as  follows: 


DPRAM 

Drag  | 

P 

J2 

J3 

J4 

J5 

fi  i 

J31 

J41 

J22 

J32 

J42 

OPRAM 

J33 

J43 

J44 

Li 

hi 

L 

X22  | 

X32 

A42 

X33 

X43 

X44 

RPRAM  -  Radar  parameters.  The  first  two  boxes  of  each  line  contain  the 
station  name.  The  succeeding  boxes  indicate  the  parameters  desired.  They 
are  ordered  as  follows; 


L 

L 

l 

|  A 

R 

A 

E 

R 

0 

This  list  applies  if  the  full  set  is  used.  See  Appendix  A  for  explanation  of 
standard  set. 
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where: 


L  =  Latitude 
i  =  Longitude 
A  =  Altitude 
R  =  Range  Bias 
A  =  Azimuth  Bias 
E  =  Elevation  Bias 
R  =  Range  Rate  Bias 
R  =  Time  Bias 


Additional  cards  may  be  added  for  more  radar  stations.  Note:  The  number 
of  X's  in  CPRAM  +  DPRAM  -f  QPRAM  must  be  <  15.  The  total  number  of 
X1  s  must  be  <  30 . 

Lines  4  1  through  46  -  Bounds 


41 

BNDS 

100 

42 

2 

100 

43 

3 

.  1 

44 

4 

45 

5 

46 

6 

A  bound  must  be  entered  for  each  parameter  selected  above  in  the  same 
sequence.  For  each  iteration  of  the  differential  correction  process,  the 
change  in  each  parameter  is: 


a.  Less  (in  absolute  value)  than  the  corresponding  bound  if 
this  bound  is  positive. 

b.  Zero  if  the  corresponding  bound  is  zero. 

c.  Unrestricted  if  the  corresponding  bound  is  negative. 

Lines  49  through  65  -  Sigmas 


49  1 

r 

SIGMA 

100 

2 

.  05 

51 

3 

.  05 

52 

12 

200 

53 

13 

.  05 

54 

14 

.  1 
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Sigmas  are  the  weighting  factors  for  the  radar  data.  Entered  here  are  a  set 
of  radar  sigmas  for  R,  A,  E,  R,  P,  Q,  P,  Q,  u,  v,  r  in  that  order,  in  feet, 
degrees,  and  seconds.  Ten  sets  may  be  entered,  1=0,  I,  2,  .  .  .  ,  9.  This 
value  of  I  is  the  one  to  be  entered  on  the  Station  Location  Card,  column  5. 
Additional  cards  may  be  inserted  here  if  necessary. 

5.  2.  3.  2  Optional  Input 

Line  29  -  Tracking  Termination  Time 

[ _ 29  H  ■  Ti  T~  ~  i  ~~ 

If  PRT IM  (21)  is  not  zero,  then  only  observations  prior  to  this  time  (in 
minutes  from  midnight)  will  be  used  in  the  orbit  determination. 


Line  5?  -  Maximum  Number  of  Iterations 


57  Ml 


MAX  IT 


If  the  differential  correction  process  has  not  converged  at  the  end  of  MAXIT 
iterations,  the  run  will  stop. 


Line  58  through  60  -  Data  Tape  Specifications 


58 

I  ’ 

IBCDI 

59 

I 

IBINI 

60 

_L 

IBINU 

_ 

IBCDI.  If  the  radar  observation  and  station  location  information  is  to  come 
in  on  a  BCD  tape  other  than  A3  (the  normal  FORTRAN  system  input  tape),  the 
tape  number  must  be  specified  here.  It  may  be  any  channel  A  tape  not  used 
by  the  system;  only  the  numeric  designation  is  required. 


IBINI.  If  TRACE  is  to  input  a  binary  tape  containing  compacted  radar  data 
(produced  by  a  previous  run),  IBINI  must  be  non-zero.  If  IBINI  <  5,  the 
tape  is  assumed  to  be  on  B5;  if  IBINI  >  5,  it  is  assumed  to  be  the  tape 
number  (Channel  B  only). 
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IBINU.  If  IBINU  is  non-zero,  TRACE  will  produce  a  binary  tape  containing 
the  sorted  processed  radar  observation  data  for  later  use.  The  same  tape 
numbering  convention  holds  as  for  IBINJ.  above.  If  IBINU  is  non-zero,  after 
a  tape  is  produced  IBINI  will  be  set  for  the  corresponding  tape  unit;  successive 
cases  of  the  same  run  will  therefore  not  require  that  the  observation  input  be 
repeated.  (The  tracking  data  input  function  must  still  be  selected  by  ITIN  in 
order  to  read  the  tape.  ) 


Line  61  -  Refractivit\ 


The  equation  used  to  correct  elevation  data  is: 


>0.  1  radian,  and  E  =  E' 

radian,  and  n  .  4  0. 

si 


1 

1000 


n  .  x  10^ 
si _ 

12  +  1000  E' 


E  =  E1  -  n  .  cotn  E1  if  E' 

si 


if  E'  <  0.  1 


E'  is  the  input  elevation.  The  n  are  read  here,  i  =  0,  1,  2,  .  .  .  ,  9.  This 

S  l 

value  of  i  is  the  one  to  be  entered  on  the  Station  Location  Card,  column  6. 

Additional  cards  may  be  inserted  here  if  necessary.  Nominally, 

n'  =  312.  0  x  10"6. 
so 


Lines  63  and  64  -  sd)S 


SOS  contains  up  to  nine  station  names  (2  blocks  per  name).  For  each  of 
these  stations,  the  root  mean  square  (rms)  and  the  square  root  of  the  sum 
of  the  squares  (SOS)  of  the  residuals  after  division  by  the  radar  sigmas  are 
printed  out. 
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Lines  65  and  following  -  Constraint  Matrix 


65 

KNST 

4 

66 

BL1ST 

1 

1 

1 

2 

2 

1 

f 

3 

)* 

3 

1 

hr 

4 

a 

4 

i 

i 

5 

l 

1 

-1 

i 

6 

i 

2 

.  5 

i 

7 

i 

5 

1 

These  quantities  can  be  best  explained  by  example.  Assume  that  there  are 
n  parameters  to  be  solved  for  (p^  p^,  .  .  .  ,  p^)  =  p.  The  ordering  of  the  p. 
corresponds  to  the  order  of  the  X's  in  CPRAM,  DPRAM,  and  RPRAM.  Also 
assume  that  there  are  m  linear  constraints  to  be  placed  on  these  parameters. 
For  example,  if  n  =  6,  m  =  2,  these  might  be  p^  +  =  6,  -  2p^  =  0.  Then 

KNST  is  equal  to  the  number  of  effective  (unconstrained)  parameters,  or 
4  { =  n-m). 


5-17 


BUST,  the  constraint  matrix,  is  obtained  as  follows: 


a.  State  the  problem  in  the  form  p  =  Bp  +  c,  where  the  p  are 
the  effective  parameters.  For  the  example  given,  this 
takes  the  form 


’pl~ 

1  0  0  0  * 

pl 

0 

P2 

0  10  0 

p2 

0 

p3 

0  0  10 

p3 

0 

• 

+ 

p4 

0  0  0  1 

p4 

0 

p5 

-10  0  0 

6 

.  p6. 

0  .5  0  0 

0 

P 


(78) 


b.  Input  the  non-zero  elements  of  the  augmented  (n  +  1)  by 
(m  +  1)  matrix 


B  c 

0  1 

where  the  element  b^j  is  input  as  i,  j,  b^j.  The  input  for 
this  example  is  shown  above. 

Lines  73  and  74  -  Time  Delay  Correction 


73 

i 

C 

74 

3 

If  C(3)  is  non-zero,  a  time  correction  will  be  applied  to  the  radar  data.  C(3) 
should  contain  ±  the  speed  of  light  (earth  radii/min). 


Then,  t' 


t  + 


R 

c[3T 


where  R  is  the  range. 


Line  75  -  Data  Editing 
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y 

X 

I ^*n - ze  ro,  the  radar  data  will  be  edited  on  the  second  and  following 
iterations.  DaTaxpptnts  will  be  discarded  with  residuals  greater  than: 

a.  The  input  sigma  times  |  C ( 1 3 )  |  if  C(13)  is  negative. 

b.  The  statistical  sigma  from  the  previous  iteration  times 
C(13)  if  3(13)  is  positive.  (A  sigma  is  calculated  for  each 
station  and  data  type.  ) 

The  above  input  should  be  followed  by  an  END  BASIC  card.  If  the  radar  data 
is  input  from  cards.  Station  Location  Data  and  Radar  Observation  sheets 
should  be  filled  out  followed  by  an  END  DATA  card.  If  the  data  is  input  from 
BCD  or  binary  tape,  the  END  DATA  card  follows  the  END  BASIC  card. 

5.  2.  3.  3 

Column 

1-? 

,  5 

6 

9-  17 

19-27 

29-36 

39-39 
(4  1-42) 


The  last  Station  Location  Card  must  be  followed  by  a  card  with  the  letters  TR 
in  columns  1-2.  There  may  be  up  to  50  stations  entered. 


Station  Location  Data 


ST.  Two  letters  that  serve  as  identification  for  a  station. 

To  t'.'r,  •■itr’.ti.r*' s  should  have  the  same  symbol 

Type  of  radar  observation  sigma  to  be  applied  to  data  from 
this  station. 

The  sets  of  sigmas  input  with  the  Basic  Data  are  numbered  (from 
0  to  9)  in  the  order  in  which  they  are  read  in.  (See  Line  49.  ) 

Type  of  refractivity  correction  to  be  used  for  elevation  readings 
from  this  station. 

R.efractivities  are  numbered  in  the  order  they  are  input  in  the 
Basic  Data.  (See  Line  61.) 

North  latitude  of  the  station  in  degrees 

East  longitude  of  the  station  in  degrees 

Altitude  of  the  station  in  feet 

If  this  station  reports  P  or  P  data  (Q  or  Q),  these  columns  contain 
the  two  letter  symbols  for  the  associated  station(s)  of  the  tracking 
configuration.  Each  such  associated  station  must  appear  on  a 
separate  Station  Location  Caid,  but  it  is  not  necessary  for 
columns  38-42  to  be  filled  out  on  the  latter. 


5-19 


5.  2.  3.4 
ST 

GMT 

DAY 

HR 

MIN 

SEC'S 

TY 

TY 

1 

2 

3 

4 

5 

6 

7 

8 
9 


Radar  Observation  Data 


Station  call  letters,  which,  must  correspond  to  a  Station  Location 
Card 

Number  of  hours  to  be  added  to  the  observation  time  to  give 
Greenwich  Mean  Time 

The  values  indicate  the  time  of  the  corresponding  observations 


Type  of  observation 

Col.  26-39 

Range 

Right  Ascension 
L 

Hour  Angle 

Af 

LI 

R 

P 

P 


Col.  41-54 
Azimuth 
De  clination 

VI 

r>e  clination 

At 

L2 

Q 

Q 


Col.  56-69 

Elevation 

N 

L3 


The  last  Radar  Observation  Card  must  be  followed  by  a  card  with  the  letters 
TR  in  columns  1-2.  (An  END  DATA  card  must  also  follow  if  no  nonstandard 
input  is  included.  ) 

5.2.  3.5  Flocking  Option 

For  large  numbers  of  radar  observations  (>1000)  the  data  should  be  divided 
into  "Flocks.  "  Flocks  may  be  of  arbitrary  size  (but  each  <1000  observations). 
A  control  card  with  the  letters  T  E  in  columns  1  and  2  is  used  to  signal  the 
end  of  a  flock;  any  number  of  these  may  be  placed  among  the  observation 
cards.  The  last  flock.  must  still  be  terminated  by  a  TR  card. 
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There  are  two  restrictions  to  be  observed.  First,  the  observations  must 
be  in  partial  chronological  order.  That  is,  every  data  time  of  a  given  flock 
must  be  later  than  all  times  in  all  previous  flocks.  Second,  the  Basic  Data 
quantity  IBINU  (Line  60)  must  be  specified  so  as  to  produce  a  compacted 
data  tape,  unless  such  a  tape,  produced  on  a  previous  run,  is  being  used 
as  input. 

The  mechanics  of  this  option  are  as  follows.  The  radar  observations  are 
read,  sorted,  processed,  and  written  on  tape,  one  flock  at  a  time,  by  TRAIN. 
If  more  than  one  flock  is  found  to  be  present,  the  differential  correction 
process  in  MAIN  reads  the  tape  a.nd  computes  residuals  and  the  normal 
matrix  for  one  flock  at  a  time. 

It  is  very  strongly  recommended  that  large  sets  of  data  be  broken  into 
flocks;  looping  and  strange  halts  result  from  overrec.ding  observations  with 
the  program. 
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Table  4.  Tracking  Input  Description 


1: 


2-8: 

9: 

10-15: 


16: 

17: 


18: 


21: 

22-28: 


29: 

33: 


Itinerary.  Specify  order  of  compu¬ 
tations. 

1-Tracking  data  input,  2-Track, 

3- Trajectory  only, 

4- Data  generation, 

5- Error  analysis. 

Epoch 

Type  1  2  3  4  of  initial  conditions 

1C  x  a  a  T 

y  6  e  6 

z  6  i  B  (ft-  sec  -deg) 
x  AO  A 
y  r  cm  r 
z  v  t  v 

NOTE:  If  r  is  negative  it  is  inter¬ 
preted  as  height;  if  v  is  neg.  ,  cir¬ 
cular  velocity  is  computed  and  used. 
CDA/W.  (ft^/lbl 
Atmosphere  Model  Specification 
0  ARDC  59 

1  Lockheed 

2  Paetzold  11 

3  Paetzold  1 

4  L.  F.  E. 

F  -  solar  radio  flux 

Ap  -  planetary  magnetic  index 

g(a)  -  plasma  intensitv  coeff. 


o  10  ° 

s!«< 


<5  2 


.H,-  -2  < 


w 

o  « 
O'  p 
s  «  <0  II  •- 
_  *-<  m  •J 


I  c  -o  ° 
X  o  (X  ii 

4)  O  3  N 


o>  a 

o 


X  to  print  £ 


II 


rm 


Print  at  t»t  (At. )t,  (At?).  .  -(^t^Jt  in 
minutes  from  midnight  if  PRT1M=1, 
from  epoch  if  0. 


23  -  n  (<9),  24  -  tQ,  25  -  Atj 
26  -  tlf  27  -  At 2,  28  -  t2,  etc. 


Fit  only  observations  prior  to  tf 
Initial  Condition  Parameter 
Specification 

Type _  X  to  solve  for: 

1  x  y  z  x  y  z  t 

2  a  6  B  A  r  v  t 

3ae  i  0  u>  T  t 

o 


I  I  I  U  HD 


34-3  5:  Differential  Eq.  Parameter 
Specification: 

X  to  solve  for: 


'■‘D  * 

~~W~  *  J2  J 3  3 4  J5  J21  J31  J41  J22 
DPRAM1  1  M  |  3  _|  j  llll 

J32  J42 


* 

J33  J43  J44  L21  L31  L41  L22  L32  L42 

opram{— )  r  t  i  r  — 1~  1  r 


'33  43  44 


rzrzi 


*or  shorter  lists  depending 
on  OBJT 

36-40:  Radar  Parameter  Specification. 

ST  L/ARAERT 


X  to  solve  for:  11  1  1  I  I  I  i  I  I H 

ST:  Station  Symbol  R: Range  Bias 
L:  Station  Latitude  A-.Azimuth  Bias 
Station  Longitude  E:  Elevation  Bias 
A:  Station  Altitude  ft:  Range  Rate  Bias 
T:  Time  Bias 

41-48:  A  bound  must  be  provided  for  every  para¬ 
meter.  For  each  iteration  of  the  differ¬ 
ential  correction  process,  the  change  in 
each  parameter  (a)  Is  (in  absolute  value) 
less  than  the  corresponding  bound  if  said 
bound  is  positive,  (b)  Is  zero  if  the 
corresponding  bound  is  zero,  (c)  Is  un¬ 
restricted  if  the  corresponding  bound  is 
negative. 

49-56:  Sigmas 

(R,  A,E,R,  P,6,  P.O.u.v.r) 

1  =  0,1,2,... 


1  is  the  sigma  type  as  referred  to  by  the 
Station  Location  Data. 

57:  Maximum  Number  of  Relations. 

58-60:  Radar  Observation  Data  Tapes 
58:  BCD  Input  Tape  if  t  A3. 

59:  Non-zero,  <  5,  if  compacted  data  is 
to  be  input  on  B5. 

60:  Non-zero,  <  5,  if  compacted  data  is 
to  be  output  on  B5. 

61:  Refractivity. 

63-64:  The  names  of  no  more  than  nine  stations. 

For  each  of  these  stations,  the  root  mean 
square  (RMS)  and  the  square  root  of  the 
sum  of  squares  (SOS)  of  the  residuals  is 
printed  out. 

6  5:  Number  of  effective  (unconstrained) 

parameters  to  be  solved  for. 

66-72:  Constraint  Matrix  Input 
74:  Time  delay  correction. 

75:  Data  editing  option. 
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Table  5.  Radar  Station  Location  Data 
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5.  2.  4 


Data  Generation 


5 .  2.  4 .  1  Required  Input 

The  only  required  input  in  addition  to  the  Basic  Data  is  the  station  location 
data  and  data  on  Data  Specification  Sheets  I  and  II.  These  will  be  described 
below. 
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5.  2.  4.  2 


Optional  Input 

Line  34  -  Rise  and  Set  Times 


1 

IFLAG 

_ 34j| X 

6 

If  IFLAG(6)  is 

0  all  data  will  be  printed  {see  sample  output  paragraph 

5. 4.  5}. 

1  rise  and  set  times  only  will  be  printed;  Data  Specification 

II  not  necessary. 

Rise  and  sets  are  at  minimum  elevation  angle  entered  on  Data  Specification 
Sheet  I. 


Line  35  -  Input  Control  for  Multiple  Cases 


35 


I  7 


Station  Location  cards  and  Data  Specification  cards  are  always  read  when  a 
4  or  a  5  is  first  encountered  in  the  1TIN  list.  In  each  following  instance  in 
the  same  I  TIN  sequence: 

if  IFLAG(7)  is 

0  Neither  Station  Location  or  Data  Specification  cards  are 

input  (same  as  previous  case). 

1  Data  Specification  is  input,  but  Station  Locations  are  not. 

-  1  Both  Station  Location  and  Data  Specification  are  input. 


Line  36  -  Observation  Tape  Generation 


36  ■  I  ETAPE  7 


If  ETAPE  is  non-zero,  then  a  BCD  radar  observation  tape  will  be  generated 
on  tape  unit  number  ETAPE.  The  tape  format  will  be  that  of  the  tracking 
input  data. 
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Lines  37  through  52  -  Noise 


37 

BB8BSI 

0  1  2  3  4  0  0  5 

38 

D 

RPRAM 

B 

0 

IDS 

s 

w 

■ 

■ 

39 

D 

■ 

■ 

Mm 

■ 

1 

■ 

l 

40 

D 

■ 

■ 

am 

■ 

■ 

■ 

■ 

41 

PBIAS 

.  001 

2 

30 

3 

-.  005 

4 

.  001 

49 

SIGMA 

50 

2 

HUB 

3 

.  01 

4 

If  NOISE  is  non-zero,  the  observations  on  the  above  tape  and  in  the  printed 
data  generation  output  will  contain  normally  distributed  random  noise  with 
mean  value  given  in  PBIAS  (Lines  41-ff)  and  standard  deviations  given  in 
SIGMA  (Lines  49-ff).  (NOISE  starts  random  number  generator.) 


If  bias  noise  is  to  be  used,  RPRAM  (Lines  38  through  40)  contain  the  radar 
observations  to  be  biased.  The  first  two  boxes  of  each  line  contain  the  station 
name.  The  succeeding  boxes  indicate  the  parameters  desired.  They  are 
ordered  as  follows: 


where: 


RPRAM 


A 


R 


R  =  Range  Bias 
A  =  Azimuth  Bias 
E  =  Elevation  Bias 
R  =  Range  Rate  Bias 
T  =  Time  Bias 


Additional  cards  may  be  added  for  more  radar  stations. 
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PBIAS  (Lines  41  -  ff.  )  contain,  the  bias  to  be  added.  The  ordering  is  the 
same  as  the  X's  in  the  RPRAM  boxes  above. 

SIGMA  (Lines  49  -  ff.  )  contain  the  standard  deviation  of  the  random  noise  to 
be  added.  Entered  here  are  a  set  of  radar  sigmas  consisting  of  values  for 
R,  A,  E,  R,  P,  Q,  P,  Q,  u,  v,  r  in  that  order,  in  feet,  degrees,  and  sec¬ 
onds.  I  sets  may  be  entered,  1  =  0,  1,  2,  .  .  .  ,  9.  This  value  of  I  is  the  one 
to  be  entered  on  the  Station  Location  Card,  column  5.  Additional  cards  may 
be  inserted  here  if  necessary. 

Line  57  -  F-ef ractivity 

1  57  B  1refr~T  1 

The  computed  elevation  is  altered  to  account  for  refraction  using  the  following 
formula: 


and 


E'  =  E  +  ns,cotn  E  if  E  >  0.  1  radian 


(79) 


E*  =  E  + 


1 


'n  .  x  10 

si 


80 


1000U2  +  1000k  *  6  +  1000E 


if  E  <  0.  1  radian  and  n  .  4  0. 

si 


(80) 


E  is  the  computed  elevation.  The  n  .  are  read  here,  i  =  0,  1,  2,  ...  9. 

This  value  of  i  is  the  one  to  be  entered  on  the  Station  Location  Card,  column 

6.  Additional  cards  may  be  inserted  here  if  necessary.  Nominally,  n  = 

-6  so 

312.  0x10  . 
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Line  65  -  Observational  Variances 


If  the  standard  deviations  in  the  six  observational  quantities  R,  A,  E,  R, 

A,  E,  are  desired,  then 

a.  The  appropriate  box  must  be  checked  on  Data  Specification 
Sheet  II, 

b.  A  covariance  matrix  for  trajectory  parameters  must  be 
supplied  in  lower  triangular  form  beginning  at  ATA(501) 

c.  The  corresponding  parameters  must  be  indicated  (as  on 
lines  37-39  of  the  Basic  Data  Input  sheet)  in  the  CPRAM, 
DPRAM,  and  OPRAM  boxes. 

The  above  input  should  be  followed  by  an  END  BASIC  card  and  an  END  DATA 
card.  Station  Location  cards  and  Data  Specification  cards  follow. 

5.  2.  4.  3  Station  Location  Data 


Column 

1-2  ST.  Two  letters,  which  serve  as  identification  for  a  station. 

No  two  stations  should  have  the  same  symbol. 

5  Type  of  radar  observation  sigma  to  be  applied  to  data  from 
this  station 

The  sets  oi  sigmas  input  with  the  Basic  Data  are  numbered 
(from  0  to  9)  in  the  order  in  which  they  are  read  in.  (See  Line 

49.) 

6  Type  of  refractivity  correction  to  be  used  for  elevation  readings 
from  this  station. 

Refractivities  are  numbered  in  the  order  they  are  input  in  the 
Basic  Data.  (See  above.) 

9-17  North  latitude  of  the  station  in  degrees 

19-27  East  longitude  of  the  station  in  degrees 
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29-36  Altitude  of  the  station  in  feet 

38-39  If  this  station  reports  P  or  P  data  (Q  or  Q) ,  these  columns 

(41-42)  contain  the  two  letter  symbols  for  the  associated  station(s) 

of  the  tracking  configuration.  Each  such  associated  station 
must  appear  on  a  separate  Station  Location  Card,  but  it  is  not 
necessary  for  columns  38-42  to  be  filled  out  on  the  latter. 

The  last  Station  Location  Card  must  be  followed  by  a  card  with  the  letters 
TR  in  columns  1  and  2.  There  may  be  up  to  50  stations  entered. 


5.  2.  4. 4  Data  Specification 

Column  Load  Sheet  I 


1-2 

9-16 

18-23 

25-30 

32-40 

5  1-58 

60-67 


Station  Call  Letters 

These  must  correspond  to  the  letters  on  some  Station  Location 
card 

Interval,  in  minutes,  at  which  data  for  this  station  is  to  be 
generated;  also  testing  interval  for  Rise-Set-only  option 

Minimum  elevation  at  which  the  vehicle  is  visible 

Maximum  elevation  at  which  the  vehicle  is  visible 
(Zero  value  will  be  set  to  90°) 

Maximum  range  (in  nautical  miles)  to  which  vehicle  is  visible 
(Zero  value  causes  this  test  to  be  ignored) 

Start  time,  from  midnight  of  start  date 
(Zero  value  implies  epoch  is  start  time) 

5  1-  52  days 
54-55  hours 
57-58  minutes 

Stop  time,  from  midnight  of  start  date 
60-61  days 
63-64  hours 
66-67  minutes 


The  last  card  of  this  type  must  be  followed  by  a  card  with  the  letters  TR  in 
columns  1  and  2. 


Column  Load  Sheet  II  (Note;  This  is  not  used  for  the  Rise-Set-only 

option) 

1-2  Station  Call  Letters 

These  must  correspond  to  the  letters  on  some  card  from  Load 
Sheet  I 
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An  X  in  the  appropriate  column  will  cause  the  quantity  listed 
above  that  column  to  be  output.  (Only  columns  7  through  14  will 
be  written  on  ETAPE  if  that  option  is  used.)  (See  Line  36.  ) 

Range  (n  mi  but  written  on  ETAPE  in  ft) 

Azimuth  (deg) 

Elevation  (deg) 

Range  Rate  (ft/sec) 

P  Dot,  Q  Dot,  P,  Q  -  Doppler  Data 

Azimuth  Rate  (deg/min) 

Elevation  Rate  (deg/min) 

Range  Acceleration  (ft/sec^) 

Mutual  Visibility 

Output  will  be  a  list  of  numbers  of  the  stations  that  are  visible 
at  the  output  time. 

(Stations  are  numbered  in  the  order  they  are  input  on  Station 
Location  cards.  )  There  is  a  maximum  of  8  stations. 

19  Latitude  of  vehicle  (deg) 

20  Longitude  of  vehicle  (deg) 

21  Surface  Range  from  station  to  subvehicle  point  (n  mi) 

22  Altitude  of  vehicle  (n  mi) 

The  following  options  require  special  input  prior  to  the  END  cards: 

23  Doppler  Rate  =  K  X  Range  Rate 
K  is  input  into  C(29) 

24  Look  Angle 

This  is  the  angle  between  an  axis  in  the  vehicle  and  the  line  of 
sight  from  the  station  to  the  vehicle.  The  direction  cosines  of 
the  vehicle  axis  must  be  in  G(37),  C(38),  and  C(39).  These  may 
be  input  as  constant,  or  the  user  may  provide  a  subroutine  called 
FANG  that  computes  the  direction  cosines  at  each  output  point. 

25  Observation  Variances 

The  standard  deviations  of  the  six  observational  quantities  R,  A, 

E,  R,  A,  E,  are  output.  These  are  based  on  a  variance- 
covariance  matrix  for  trajectory  parameters.  The  matrix  is 
input  beginning  in  ATA(50  1)  in  lower  triangular  form  (see  Line  69). 
Corresponding  parameters  must  be  indicated  in  the  CPRAM, 
DPRAM,  and  OPRAM  boxes  (see  Line  65). 


n 

\ 


7-25 

7 

8 

n 

/ 

10 

11-14 

15 

16 
’.7 
18 
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Table  7.  TRACE  -  Data  Generation  Input 


X-1  7090  INPUT  DATA 


AEROSPACE  CORPORATION 

COMPUTATION  &  DATA  PROCESSING  CENTER 


Table  8. 


Data  Generation  Input  Description 


1:  Itinerary.  Specify  order  of  compu¬ 

tations. 

1-Tracking  data  input,  2-Track, 

3- Trajectory  only, 

4- Data  generation, 

5- Error  analysis 

2-8:  Epoch 

9".  Type  1  2  3  4  of  initial  contitions 

10-15:  IC  x  o  a  T 

y  6  e  6 

z  P  i  8  (ft-sec-deg) 

^  A  0  A 
y  r  o:  r 

z  v  t  v 


NOTE:  If  r  is  negative  it  is  inter¬ 
preted  as  height;  if  v  is  neg.  ,  cir¬ 
cular  velocity  is  computed  and  used. 
16:  Cj^A/W  (ft2/lb) 

17:  Atmosphere  Model  Specification 


18: 

19: 

20: 


0  ARDC  59 

1  Lockneed 

2  Paetzold  II 

3  Paetzold  1 

4  L.  F.  E. 

-  solar  radio  flux 

-  planetary  magnetic  index 

-  plasma  intensity  coeff. 


>•  3  in 

u  J2  u  .  w  *  o  u 

2  2  ^  c  c  o  £ 

.Z  —  <J  .  r  w  rt  ••  • 

'  u  t  5  C-OO  » 

h  -  o  a  II  ^ 

O  3  H  O 


<u  rt  .  . 

I-  Q-<  >  O 


21:  X  to  print  L1..1.I  LI  DTUD 


22-28:  Print  at  t  =  t  (At.  )t.  (At-, ).  .  .  (At  )t  in 
o  1  I  (.  n  n 

minutes  from  midnight  if  PRTIM-1, 
from  epoch  if  0. 

23  -  n  (  <9),  24  -  t  ,  25  -  At. 

—  o  i 

26  -  tj ,  27  -  Atz,  28  -  t2>  etc. 

34:  If  non-zero,  rise  and  set  times  only 

will  be  generated. 

35:  Input  control  for  station  location 

and  ephemeris  cards. 

36:  Tape  number  of  BCD  radar  station 

and  observation  tape.  (If  zero,  no 
tape  is  generated.  ) 


37:  If  non-zero,  the  above  tape  will 

contain  normally  distributed  random 
noise  (mean  values  in  PBIAS  and 
standard  deviations  in  SIGMA;  RPRAM 
specifies  the  bias  parameters). 
38-40:Radar  Parameter  Specification 

C  T>  T>  A  TP  i>  T 

X  to  specify  |  1  I  Q)  ~L~  1  1  1  1  1 

ST:  Station  Symbol  R:  Range  Bias 

A:  Azimuth  Bias 
E:  Elevation  Bias 
R:  Range  Rate  Bias 
T:  Time  Bias 

41 -48:Contains  the  mean  values  of  the  biases 

to  be  added  to  the  observations  specified 
in  RPRAM. 

49-56:Contains  the  standard  deviations  of  the 
random  noise  to  be  added. 

57-64;  Ref ractivity. 

o5:  Initial  Condition  Parameter  Specification: 

TYPE  [X]  to  specify: 

1  X  y  z  x  y  z  tQ 

2  a  6  (3  A  r  v  tQ 

3  a  e  i  ft  o  t  tQ 


66-67-Differential  Eq.  Parameter  Specification: 

X  to  specify: 


D  •  * 

1  "W  1  V-  J2  J3  J4  J5  J2 1  J31 


DPRAMl  |!l  1  I 

.  1  L  L 

J41  J22  J 

32  J42 

1  1  ri 

J33  J43  J44  T<2 1 

E31 

L41 

opRamI _ 1  1  1  1  1  I 

T*22  L32  L42 

L33 

L43  L44 

T  T  T  1  f  1  1 

6 9-76:Covariance  matrix  in  lower  triangular 
form  for  parameters  selected  above. 

••^This  list  applies  if  <t>BJT  =  4.  If  <j)BJT  =  3 
replace  with  shorter  list  J31,  ^22’  ^32'^33 

and  j  ,  L31,  E32,  E33.  If  (t)BJT  =  2 


replace  with  J^2  and  E^j,  E 2Z • 
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Table  9-  Radar  Station  Log 


Table  8.  Data  Generation  Input  Description 


Itinerary.  Specify  order  of  compu¬ 
tations. 

1-Tracking  data  input,  2-Track. 

3- Trajectory  only, 

4- Data  generation, 

5- Error  analysis 

Epoch 

Type  1  2  3  4  of  initial  contitions 


1  0-  1  5:  IC  x  ct  a  t 
y  6  e  6 

z  S  i  0  (ft-sec-deg) 
x  A  C.  A 
y  r  u:  r 
z  v  t  v 

NOTE:  If  r  is  negative  it  is  inter¬ 
preted  as  height;  if  v  is  neg.  ,  cir¬ 
cular  velocity  is  computed  and  used. 
16:  CjyV/W  (ft2 /lb) 

17:  Atmosphere  Model  Specification 

0  ARDC  59 

1  Lockheed 

2  Paetzold  11 

3  Paetzold  I 

4  L.  F.  E. 

18:  F  -  solar  radio  flux 

19:  A^  -  planetary  magnetic  index 

20:  g{a)  -  plasma  intensity  coeff. 

H 

u'  jo  o  .  0)  2  o  ” 

2  2,  2  o-  c  c  ^  £ 

«  2  s  <  ~ 

«  “  i  coo  “ 

Z  u  a<  >ou3^o 

21:  X  to  print  (  I  I  I  I  I  I  I  \  T  I  I  1 

22-28:  Print  at  t-t  (At  )t,  (At? ).  .  .  (At  )t  in 
o  1  1  n  n 

minutes  from  midnight  if  PRTIM=1, 
from  epoch  if  0. 

23  -  n  {  <9),  24  -  t  ,  25  -  At, 


37:  If  non- zero,  the  above  tape  will 

contain  normally  distributed  random 
noise  (mean  values  in  PBIAS  and 
standard  deviations  in  SIGMA;  RPRAM 
specifies  the  bias  parameters). 
38-40:Radar  Parameter  Specification 

QT  RATT'RT  - 

X  to  specify  1  I  1  ©  I  1  I  I  II 

ST:  Station  Symbol  R:  Range  Bias 

,  ,  A:  Azimuth  Bias 

©blank  spaces  £.  Elevation  Bias 

R:  Range  Rate  Bias 
T:  Time  Bias 

41  -48:Contains  the  mean  values  of  the  biases 

to  be  added  to  the  observations  specified 
in  RPRAM. 

49- 56: Contains  the  standard  deviations  of  the 
random  noise  to  be  added. 

57-64:  Refractivity . 

o5:  Initial  Condition  Parameter  Specification: 

TYPE  [X]  to  specify: 

1  x  y  z  x  y  z  tQ 

2  a  6  (3  A  r  v  tc 

3  a  e  i  ft  u  r  t0 


66-b7:Differential  Eq.  Parameter  Specification: 

X  to  specify: 


p  J2  J3  J4  J5  J2i  J31 


DPRAM 


J41  J22  J32  J42 


J33  J43  J44  L21  L31  L41 


OPRAM 


26  -  tj. 


27  -  At^,  28  -  t^.etc. 


L22  L32  l42  l33  L43  L44 


If  non-zero,  rise  and  set  times  only 
will  be  generated. 

Input  conti ol  for  station  location 
and  ephemeris  cards. 

Tape  number  of  BCD  radar  station 
and  observation  tape.  (If  zero,  no 
tape  is  generated.  ) 


6 9- 76Govariance  matrix  in  lower  triangular 
form  for  parameters  selected  above. 

*This  list  applies  if  (J)BJT  =  4.  If  0BJT  =  3 
replace  with  shorter  list  J2i>  J31,  ^ZZ‘  J32>^33 

and  L2  j  ,  L31,  L22i  L32,  ^33-  ^  tf)BJT  =  2 
replace  with  J2j,  J22  ar,d  ^21’  ■^'22* 
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Da 


•  CPRAM  -  Initial  condition  parameters.  The  first  box 
specifies  the  tvxte  of  initial  condition.  The  succeeding 
boxes  indicate  the  particular  parameter  desired.  The 
boxes  are  ordered  as  follows: 


Type 


CPRAM  1 

X 

y 

z 

X 

y 

z 

2 

a 

6 

0 

A 

r 

V 

3 

a 

e 

i 

0 

U) 

T 

DPRAM  and  (|)PRAM  -  Differential  equation  parameters. 
The  boxes  are  ordered  as  follows: 


DPRAM 

Drag  j  p 

:  J2 

J3 

J4 

J5 

* 

J21 

J31 

J41 

J  22 

J32 

J42 

<f)PRAM 

i 

I  J33  iJ43. 

l  44 

❖ 

*21 

X31 

X41 

X22 

'32 

X42 

X  3  3 

X43 

X44 

•  RPR  AM  -  Radar  parameters.  The  first  two  boxes  of  each 
line  contain  the  station  name.  The  succeeding  boxes  indi¬ 
cate  the  particular  parameter  desired.  They  are  ordered 
as  follows: 


RPR  AM 


where 


L  =  Latitude 
f  =  Longitude 
A  =  Altitude 
R  =  Range  Bias 
A  =  Azimuth  Bias 
E  =  Elevation  Bias 
R  =  Range  Rate  Bias 
T  =  Time  Bias 
u  =  ’’u"  Bias 
v  =  "v"  Bias 

This  list  applies  only  if  the  full  set  is  used.  See  the  appendix  for 
explanation  of  standard  set. 
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Additional  cards  may  be  added  for  more  radar  stations.  Note:  The  number 
P's  plus  the  number  of  Q's  in  CPRAM  +  DPRAM  -f  (j)PRAM  must  be  <15.  The 
total  number  of  P's  plus  Q's  must  be  <  30. 

Lines  37  through  44  -  Sigmas 


37 

SIGiVA 

100 

38 

2 

.  05 

39 

3 

.  05 

40 

10 

50. 

41 

12 

1  20 

42 

13 

.  Co 

43 

■ 

14 

.  06 

44 

21 

60. 

Sigmas  are  the  weighting  factors  for  the  radar  observation  partials.  Each 
sigma  is  a  standard  deviation  for  the  particular  observation  type  and  station. 

A  set  of  radar  sigmas  consists  of  values  for  R,  A,  E,  R,  P,  Q,  P,  Q, 

u,  v#  r.  in  that  order,  in  feet,  degrees,  and  seconds,  I  sets  may  be  entered, 

1  =  0,  1,  2 . 9.  This  value  of  I  is  the  one  to  be  entered  on  the  Station 

Location  Card,  column  5.  Additional  cards  may  be  inserted  here  if  necessary 


5.2.  5.  2 


Optional  Input 

Line  45  -  Covariance  Output  Specifications 

(If  any  covariance  matrix  output  is  desired  the  eighth  box,  labeled 
"update,"  of  PRCDE  must  be  checked.  ) 


45  BpIpRcibv  1x1  Id]  Ixl  I  |p|  |d|  ~lx 


An  X  in  a  box  specifies  that  the  whole  covariance  matrix  be  output.  A  D  in 
a  box  causes  only  the  square  roots  of  the  diagonals  to  be  output.  The  boxes 
are  ordered  as  follows: 
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O  <D 

0  0 

OB  Ofl 

.K. 

o 

O.  u 

-k. 

o 

c 

rt 

^  Q.  # 

c 

* 

r-H 

4-> 

a 

■  *H 

C/l 

& 

_ _ _ 

c 

“  aj 

■  rH 

w 

a. 

c 

<D 

u 

0  u 

0 

u 

4-> 

u 

rQ 

a 

r—i 

£ 

c  s 

4-> 

U 

rH 

£ 

oS 

O 

d 

o 

U 

o 

0, 

W 

&  £ 

U 

o 

& 

w 

' — ' 

' — '  ' 

■ — ■ 

**n. 

PRCOV 

0 

U 

r'-H 

£ 

u 

r-H 

N 

u  ; 

U  . 

w 

u 

p  1 
u  ; 

CM 

u 

fM 

X 

u 

CM 

hi 

U 

CM 

U 

CM 

w 

u 

i 

i 

Line  46  -  Addii 

tional  Option  Box 

C 

m  raEu5i-io>;4BlE3DI 

An  explanation  of  the  purpose  of  each  box  follows: 


$PB(pX 

L 

A. 

If 

"A1 

B. 

If 

"B1 

C. 

If 

"C1 

A 


B 


If  "A"  box  contains  an  X  the  dP/dQ  will  be  printed. 

T 

If  "B"  box  contains  an  X  the  A  A  print  will  be  omitted. 

box  contains: 

T 

1  ,  the  A  A  will  be  punched. 

T 

2,  the  partitioned  Ap  Ap  matrix  will  be  punched. 
Lines  49  through  60  -  C(Q)  Covariance  Matrix  Input 


49 

•KWf«U 

.01 

50 

2 

.  01 

51 

3 

.  01 

52 

4 

5. 

53 

5 

5. 

54 

6 

2500. 

55 

56 

These  covariance  matrices  include  the  "Q"  effects. 
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C(T)2>  (Period,  Apogee, 


COVQ  contains  the  variance-covariance  matrix  of  the  "O"  parameters,  which 
are  specified  in  CPRAM,  DPRAM,  OPRAM.  and  RPRAM.  The  matrix  is 
input  in  lower  triangular  form.  For  example,  if  there  were  3  O's  specified, 
the  matrix  would  have  the  form: 


This  matrix  is  input  into  COVO  in  the  order: 


If  IFLAG(7)  is: 

-  1  Input  all  Station  Location  and  Data  Specification  Cards 

0  No  input  -  values  are  the  same  as  previous  case 

1  Input  Data  Specification  cards  only;  station  locations  are  the 
same  as  previous  case. 

This  option  applies  to  "multiple”  cases  (that  is,  ITIN  =  553  .  .  .  ),  but  not  to 
"stacked"  cases  (successive  cases  for  each  of  which  ITIN  =  5). 


If  IFLAG  (15)  is: 


T 

-  1  Hold  A  A  from  previous  case 
T 

0  No  A  A  or  inverse  input 

1  Input  A^A  into  ATA  area,  as  augmented  upper  triangular  matrix 

2  Input  inverse  into  ATA  (501)  area,  as  lower  triangular  matrix. 

5.  2.  5.  3  Station  Location  Data 
Column 

1-2  ST 

Two  letters,  which  serve  as  identification  for  a  station. 

No  two  stations  should  have  the  same  symbol 

5  Type  of  radar  observation  sigma  to  be  applied  to  data  from  this 
station 

The  sets  of  sigmas  input  with  the  Basic  Data  are  numbered 
(from  0  to  9)  in  the  order  in  which  they  are  read  in  (see  Lines  37 
through  44). 

6  Not  used  for  Error  Analysis 

9-17  North  latitude  of  the  station  in  degrees 

19-27  East  longitude  of  the  station  in  degrees 

29-36  Altitude  of  the  station  in  feet 

38-39  If  this  station  reports  P  or  P  data  (Q  or  Q),  these  columns 

(4  1-42)  contain  the  two  letter  symbols  for  the  associated  station(s) 

of  the  tracking  configuration.  Each  such  associated  station 
must  appear  on  a  separate  Station  Location  card,  but  it  is  not 
necessary  for  columns  38  through  42  to  be  filled  out  on  the 
latte  r. 

The  last  Station  Location  card  must  be  followed  by  a  card  with  the  letters 
TR  in  columns  1  and  2.  There  may  be  up  to  50  stations  entered. 
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5.  2.  5.4 


Data  Specification 


Load  Sheet  I 


Column 

1-2  Station  Call  Letters 

These  must  correspond  to  the  letters  on  some  Station  Location 
card 

9-16  Interval,  in  minutes,  at  which  data  for  this  station  is  to  be 

gene  rated 

18-23  Minimum  elevation  at  which  the  vehicle  is  visible 

25-30  Maximum  elevation  at  which  the  vehicle  is  visible 

(Zero  value  will  be  set  to  90°) 

32-40  Maximum  range  (in  n  mi)  to  which  vehicle  is  visible 

(Zero  value  causes  this  test  to  be  ignored) 

5  1-58  Start  time,  from  midnight  of  start  date 

(Zero  value  implies  epoch  is  start  time) 

51-52  days 
54-55  hours 
57-58  minutes 

60-67  Stop  time,  from  midnight 

60-6  1  days 
63-64  hours 
66-67  minutes 

The  last  card  of  this  type  must  be  followed  by  a  card  with  the  letters  TR  in 
columns  1  and  2. 

Load  Sheet  II 


Column 

1-2  Station  Call  Letters 

These  must  correspond  to  the  letters  on  some  card  from 
Load  Sheet  I. 

7-18  An  X  in  the  appropriate  column  will  cause  the  quantity  listed 

above  that  column  to  be  computed  and  used  internally. 

The  last  card  of  this  type  must  be  followed  by  a  card  with  the  letters  TR  in 
columns  1  and  2. 


The  only  limit  on  the  number  of  cards  using  these  formats  is  that  at  most  fifty 
different  stations  are  allowed. 
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Table  13.  Error  Analysis  Input  Description 


1: 


2-8: 

9: 

10-15: 


Itinerary.  Specify  order  of  computa-  3^. 
tions. 

1 -Tracking  data  input,  2- Track, 

3- Trajectory  only 

4- Data  generation 

5- Error  analysis 
Epoch 

Type  1  2  3  4  of  initial  conditions. 

IC  x 


Extension  of  Diff.  Eq.  Parameter 
Specification 

jpj  Parameter  to  be  solved  for. 
jo|  Pararnejer  considered  in  error 

J33J43J44L21L31L41L22L32L42L33L43L44 

i.  r  l  i  1  1  .1  .  }  rnzi 


1 

6 

8  (ft-sec-deg) 

A 

r 


16: 

17: 


18: 


y  C 

2  f 

X  1 

y  * 

z  \ 

NOTE:  If  r  is  negative  it  is  inter¬ 
preted  as  height;  if  v  is  neg.  ,  cir¬ 
cular  velocity  is  computed  and  used. 
CdA/W.  (ft  2/ lb) 

Atmosphere  model  specification 
0  ARDC  59 

1  Lockheed 

2  Paetzold  II 

3  Paetzold  I 

4  L.  F.  E. 

F  -  solar  radio  flux 

Ap  -  planetary  magnetic  index 

g(a)  -  plasma  intensity  coeff. 


21.; 

22-28: 


29: 


30: 


\ 


U  » 

3  2 


to  tl  ^ 
<U  «S 

u 


—  C  ( 
u  V  « 

•  P  m  nj 
£  V  c  *0  ' 

$  -S  o  a 
>  v  u  3 


32-36:  Radar  Parameter  Specification. 

13  to  solve  for 
£51  parameter  in  error 
S  T  L  £  ARAERTuv 

LITTL  .1  ITT  i  m 

ST:  Station  Symbols  R:  Range  Bias 
L:  Station  Latitude  A:  Azimuth  Bias 
l:  Station  Longitude  E:  Elevation  Bias 
A;  Station  Altitude  R:  Range  Rate  Bias 
T:  Time  Bias 
u:  "u"  Bias 
v:  "v"  Bias 

37-44:  Sigmas 

(R,  A,  E,  R,  P,  6,  P,  O,  u,  v,  r) 

Is  0,  1 ,  2,  ....  9  1 

I  is  the  sigma  type  as  referred  to  by  the 
Station  Location  Data. 

45:  Covariance  Matrix  Output  Specification. 

155  in  box  for  whole  matrix 
(Sj  in  box  for  sq.  rt.  of  diagonals 


46: 


0 

ft 

0 

0 

hi 

U 

0 

u 

H 

O 

oj 

ft 

O 

O 

uT 

O 

U 

(M 

0 

oO 

-4-1 

H 

b 

1 

,etc. 


X  to  print  |  |  |  1  |  |  |  |  |  |  1  1  | 

Print  at  t=tQ(Atj )t^ (At^).  .  .(At  )t  in 
minutes  fro'm  midnight  if  PR^’I^=1, 
from  epoch  if  0. 

23  -  n{<  9),  24  -  t  ,  25  -  At 

26  _  tj,  27  -  A°t2,  28  -  t2 

Initial  Condition  Parameter  Specifi¬ 
cation 

[H  Parameter  to  be  solved  for. 

O  Parameter  considered  in  error. 
Type 

1  xyzxyzt 

2  06  8  A  r  v  t° 

3  aei  0  <»  t  t 

i  li.  11L  r  1 id 

Differentia]  Equation  Parameter 
Specification 

Parameter  to  be  solved  for. 
Parameter  considered  in  error. 


Option  Boxes:  The  form  is: 

OPBOX  IaIb lei 

A.  If  MA"  box  contains  an  X  the  5P/SQ 

will  be  printed.  „ 

B.  If  "B"  box  contains  an  X  the  A  A 

will  not  be  printed.  _ 

C.  If  "C"  box  contains:  1,  the  A  A  will 
be  punched;  2,  the  partitioned  A  xAp 
will  be  -punched 

49-60:  C(O)  covariance  matrix  in  lower  tri¬ 
angular  form. 

62:  Input  control  for  station  location  and 

ephemeris  cards. 

T  T  - 1 

63:  A  A  or  (A  A)  input  option. 

-1,  hold  A1  A  from  previous  case 

0,  no  ATA .or  inverse  input 

1,  input  A-^A  into  ATA  area 

2,  input  inverse  into  ATA (501)  area. 


*  or  shorter  list  depending  on  <f)BJT 


yiJ2.T3J4J5J2lJ3lJ41  J22  ^32^42 
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Table  14.  Radar  Station  Location  Data 


► 
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5.  2.  6 


Deck  Setup 


TRACE  may  be  run  from  either  binary  cards  or  tape.  This  section  describes 
the  setup,  from  the  user's  point  of  view,  of  each  mode,  indicates  how  to  pro¬ 
duce  a  program  tape,  and  explains  the  use  of  the  dummy  routines  provided 
with  the  deck. 

5.  2.  6.  1  Running  from  Tape 

At  this  point  it  is  necessary  to  introduce  REIN,  a  single  program  with  only 
one  function:  to  read  CHAIN  {the  first  link  of  TRACE)  from  some  specified 
tape  and  thus  initiate  execution  from  this  tape. 

When  a  binary  program  tape  is  available,  the  input  setup  is  as  shown  in 
Figure  16.  REIN  can  be  considered  a  loader,  which  calls  CHAIN  from 
logical  unit  8  (A-8  with  the  present  Aerospace  Unit  Table).  REIN  must  be 
reassembled  if  any  other  unit  is  desired. 

All  comments  below  pertaining  to  the  deck  organization  when  running  from 
cards  are  applicable  to  tape  also,  as  every  execution  from  cards  first  pro¬ 
duces  {and  then  uses)  a  program  tape. 

When  using  a  previously  written  tape,  it  is  not  possible  to  compile  and  then 
use  any  program  other  than  REIN. 


INPUT 


['•'DATA 
/ REIN 


7/ 

) 


Figure  16.  Running  from  Tape 
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5.  2.  6.  2  Running  from  Cards 
5.  2.  6.  2.1  Arrangement  of  Deck 

The  complete  program  consists  of  five  links:  CHAIN,  TRAIN,  MAIN,  GAIN, 
and  FEIGN.  (Use  of  REIN  applies  to  execution  from  tape  only.  )  A  standard 
all-purpose  deck  would  look  like  Figure  17. 

Arrows  indicate  positions  of  symbolic 
programs  and/or  Debug  cards  to  be 
used  with: 


INPUT  '/ 


♦DATA 


/  FEIGN 


♦CHAIN  (5,  B 3) 


GAIN 


♦CHAIN  (4,  B 3) 


MAIN 


♦CHAIN  (3,  B3)'| 


TRAIN 


7 


I  *  CHAIN  (2.  B  3) 

- 7; 

CHAIN  / 


♦CHAIN  (1 ,  B  3) 


/ 


V 


y 


FEIGN 


GAIN 


MAIN 


TRAIN- 


CHAIN 


♦XEO 


For  various  reasons,  a  varie.ty  of  smaller  decks  or  tapes  may  prove  more 
desirable.  It  is  necessary  to  include,  for  any  given  run,  only  those  links 
to  be  executed  during  that  run.  The  types  of  execution  and  the  links  that 
each  require  are: 


T  racking 
Trajectory  only 
Data  generation 
Error  analysis 


CHAIN,  TRAIN,  MAIN 
CHAIN,  MAIN 
CHAIN,  GAIN 
CHAIN,  FEIGN 


Some  economy  in  tape-handling  results  from  the  use  of  shorter  decks  in 
production  work. 


The  input  quantity,  PTAPE  =11,  must  be  included  in  the  Basic  Data  when 
running  from  cards  (see  paragraph  5.  2.  6.  2.  2). 

If  symbolic  cards  for  compilation  are  to  be  included,  they  must  immediately 
follow  the  CHAIN  control  cards  and  precede  the  DEBUG  cards  (if  any),  for 
the  appropriate  link.  In  this  connection,  it  is  usually  worthwhile,  but  not 
necessary,  to  remove  the  corresponding  binary  program  from  the  link. 

5.  2.  6.  2.  2  Producing  a  Tape 

While  it  is  true  that  every  run  from  cards  automatically  produces  a  program 
tape,  a  short  explanation  of  the  CHAIN  control  card  preceding  each  link  may 
clarify  the  process. 

The  card 

*  CHAIN  (I,  B3)  (1=1,  2,  3,  .  .  .  ) 

assigns  the  number  I  to  the  link  it  precedes,  and  directs  the  FORTRAN 
monitor  to  store  this  program  on  tape  B3. 

Execution  begins,  after  each  program  has  been  stored  on  the  tape  designated 
by  its  control  card,  by  reading  back  in  the  (physical)  first  link  of  the  deck; 
this  must,  therefore,  always  be  CHAIN.  If  the  tape  number  (in  this  case,  B3) 
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is  the  same  on  each  control  card,  a  binary  tape  of  TRACE  is  then  available 
by  simply  saving  B3.  When  source  programs  are  included  with  the  binary 
cards,  the  result  of  the  compilation  appears  on  the  tape. 


Links  are  brought  into  core  by  programming  in  TRACE,  which  assumes  them 
to  be  on  the  tape  designated  by  the  input  quantity  PTAPE.  PTAPE  is  set,  in 
REIN,  to  8,  and  does  not  have  to  be  read  during  a  normal  run  from  tape. 

When  using  cards,  however,  PTAPE  must  be  input.  It  must  equal  11  unless 
the  control  cards  are  changed  (B2  and  A4  are  the  only  other  units  FORTRAN 
will  recognize  for  this  purpose  at  present);  or  a  unit  table  is  employed  in 
which  B3  does  not  equal  11. 

5.Z.6.2.3  Dummy  Routines 

For  some  purposes,  it  may  be  desirable  to  increase  the  amount  of  core 
storage  available  for  data  handling. 

Since  TRACE  contains  many  options,  not  all  of  which  are  usually  executed 
in  any  one  run,  there  are  always  a  certain  number  of  extraneous  subroutines 
present.  (For  instance,  only  one  integration  routine  is  ever  in  use  during 
any  trajectory.)  These  routines  may  be  replaced  with  one-word  ''dummies" 
if  more  storage  cells  are  needed.  (In  replacing  cards,  it  should  be  noted  that 
the  TRACE  programs,  and  then  the  library  routines,  are  alphabetical  within 
each  link.  ) 

A  list  of  a  few  subroutines  of  appreciable  size  that  might  be  dummied  are 
given  here  as  an  example. 
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Link 

Routine 

When  It  Is  Unnecessary 

Approx.  No.  of  Words 

MAIN 

SDUMP 

No  dump  required 

150 

DXDA 

No  analytic  partials 

400 

DXDRSR 

DXDA  being  used 

350 

PATTY 

No  accumulated  normal 
matrix  output 

400 

PTRAJ 

No  trajectory  output 

170 

REST 

BLIST  (  constraints)  not  used 

260 

resti 

i<  ii  1 1  m 

215 

AMRK 

Gauss- Jackson  integration 
used  (  COW) 

390 

COW 

AMRK  used 

1  100 

5.  2.  6.  3  Data  Deck  Setup 

All  types  of  runs  require  input  of  BASIC  DATA  followed  by  the  END  BASIC 
card.  ITIN  determines  what  further  input  is  required. 


If  ITIN  = 

1 

4 

5 

(TRAIN) 

(GAIN) 

(FEIGN) 

Read 

Tracking  Data  Input 

Data  Generation 

Error  Analysis 

Station  Location 

1  St 

2nd 

2nd 

Radar  Observation 

2nd 

Data  Specification 

I  and  II 

3rd 

3rd 

END  DATA  Card 

1  St 

1  St 

ITIN  =  2,  3,  requires  only  an  END  DATA  card. 

There  is  one  exception  to  the  above  chart.  If  two  or  more  GAIN  or  FEIGN 
runs  are  run  in  the  same  ITIN  sequence,  the  Station  Location  and  Data 
Specification  cards  are  normally  read  the  first  time  only.  (See  pages  5-28 
and  5-43. ) 
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For  two  successive  trajectories  the  deck  in  Figure  18  should  be  followed  by 
the  one  in  Figure  19. 


Figure  19.  ITIN  =  33,  Two  Trajectories 


The  input  deck  for  a  standard  tracking  run  with  observation  cards  is  shown 
in  Figure  20. 


Figure  20.  ITIN  =  12,  Tracking  Input  and  Run 
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For  a  tracking  run  followed  by  data  generation,  the  above  deck  should  be 
followed  by  Figure  21. 


Figure  21.  ITIN  =  124,  Tracking  Run  Plus  Data  Generation 

Data  generations  and  error  analyses  are  run  from  input  decks  of  identical 
structure,  as  in  Figure  22. 


•j  Figure  22.  ITIN  =  4  or  5,  Data  Generation  or  Error  Analysis 
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OUTPUT 


5.  3 

Each  link,  or  type  of  computation,  chosen  in  TRACE  provides  two  types  of 
programmed  output.  First,  there  are  those  headers  and  quantities  that  are 
a  function  only  of  the  link  being  executed  and  are  not  controlled  by  input 
options.  Second,  there  is  output  which  must  be  specifically  selected  through 
the  input  quantities  PRCDE,  PRTIM,  and  some  of  the  tape  parameters.  A 
third  possible  type  of  output  is  available,  to  anyone  sufficiently  familiar  with 
the  program,  by  use  of  the  FORTRAN  DEBUG  capabilities. 

The  first  two  kinds  of  output  can  best  be  described  by  means  of  examples. 
Fifteen  sample  printouts  are  included  at  the  end  of  this  section. 

5.  3.  1  Examples  1  and  2  (CHAIN) 

The  output  from  CHAIN  is  all  of  the  first  type  and  appears  on  every  run. 

The  first  26  lines  are  the  BCD  card  images  of  the  FINP  input.  Inclusion  of 
a  card  with  the  symbol  CLOCK  in  columns  7  through  11  will  cause  the  time 
(in  hundredths  of  a  minute),  at  which  this  card  is  read,  to  be  output.  Printing 
of  the  input  cards  may  be  eliminated  by  removing  the  special  versions  of  sub¬ 
routine  (CSH)S  from  the  binary  deck.  Any  error  printout  indicates  that  the 
last  card  was  either  punched  incorrectly,  includes  a  symbolic  location  not  in 
the  FINP  list,  or  is  not  a  Basic  Data  card. 

5.  3.  2  Example  3  (TRAIN) 

TRAIN  output  is  not  input-controlled  and  is  produced  whenever  a  tracking 
run  is  executed.  Example  3  is  the  result  of  specifying  station  location  and 
radar  observations  on  cards,  along  with  the  normal  input  deck.  If  a  binary 
radar  data  tape  had  been  read,  TRAIN  output  would  not  include  the  observa¬ 
tional  data. 

There  are  several  possible  error  messages.  In  Example  3,  a  mispunched 
card,  redundancy  in  reading  the  BCD  tape,  and  inclusion  of  an  observation 
from  a  station  for  which  there  is  no  station  location  card,  would  all  produce 


the  same  effect.  An  appropriate  line  of  output  is  printed,  the  observation 
(or  station)  in  question  is  deleted,  and  execution  continued.  Redundancies 
in  reading  the  binary  data  tape  would  cause  this  information  to  be  printed 
and  execution  terminated. 

5.3.3  Examples  4,  5,  6,  and  7  (MAIN  -  Trajectory  Only) 

All  output  from  MAIN  during  a  trajectory  run,  with  the  exception  of  the 
initial  conditions,  must  be  selected  through  the  input  quantities  PRCDE  and 
PRTIM. 

Example  4  results  from  the  request  that  the  constants  in  use  in  the  program 
be  printed  out. 

Example  5  shows  trajectory  output;  Example  6,  the  variational  equations; 
and  Example  7,  the  elements.  In  these  examples,  the  three  kinds  of  output 
are  shown  singly;  any  two,  or  all  three,  could  have  been  given  at  the  same 
time  just  as  easily.  However,  it  is  not  possible  to  request  one  type  of  output 
at  one  sequence  of  times  and  another  type  at  another  sequence  during  the 
same  trajectory. 

5.3.4  Examples  8,  9,  10,  and  11  (MAIN  -  Tracking) 

When  using  TRACE  as  a  tracking  program,  MAIN  produces  both  input- 
independent  and  input-controlled  output. 

The  first  type  consists  of  the  initial  conditions,  and,  for  each  iteration, 
something  similar  to  Example  8.  The  legend,  CURRENT  SOLUTION  IS  NOT 
GOOD,  indicates  that  the  previous  solution  has  caused  the  rms  of  the  residuals 
to  increase.  The  program  therefore  will  decrease  the  bounds,  return  to  the 
last  good  solution,  and  re-solve,  using  the  corresponding  normal  matrix. 
CURRENT  SOLUTION  IS  BEST  SO  FAR  is  a  signal  that  the  rms  has  decreased 
and  that  the  bounds  will  be  increased  for  faster  convergence.  SIGMA 
(PARAMETERS) /SIGMA  (NORMALIZED  DATA)  is  the  square  root  of  the 
diagonal  of  the  inverse  normal  matrix. 
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The  input-controlled  output  includes  everything  that  can  be  obtained  during 
a  trajectory-only  run  (see  paragraph  5.  3.  3),  plus  four  additional  computa¬ 
tions  . 

Example  9  shows  the  STATION-BY-STATION  SOS  (square  root  of  the  sum 
of  the  squares  of  the  residuals  after  division  by  the  appropriate  radar  sigmas). 
This  is  printed  once  per  iteration  and  is  limited,  at  present,  to  nine  stations. 

Example  10  gives  the  measured-minus-computed  values  of  the  radar  residuals. 
If  a  time  bias  is  included  from  one  or  more  stations,  the  output  time  will  be 
biased. 

The  partials  of  the  radar  observations,  with  respect  to  the  parameters  being 
solved  for,  appear  as  in  Example  11.  This  may  be  requested  along  with,  or 
independently  of,  the  residuals  (Example  10).  The  units  of  the  partials  are 
in  earth-radii  and  radians;  (all  other  output  is  in  feet  and  degrees,  or  any 
other  system  specified  using  the  non-standard  input);  these  quantities  are 
output  before  the  division  by  the  radar  sigmas. 

5.  3.  5  Examples  12,  13,  14,  and  15  (GAIN  -  Data  Generation) 

All  the  output  illustrated  in  paragraph  5.  3.  3  is  also  available  during  a  data 
generation.  Besides  this,  there  are  three  additional  types  of  output; 

Example  12  shows  the  station  locations  and  data  specifications. 

Example  13  is  a  result  of  choosing  to  employ  GAIN  in  the  Ri se -Set-Only 
mode . 

Example  14  shows  the  various  types  of  data  that  the  program  can  produce. 

5.  3.  6  Example  15  (FEIGN  -  Error  Analysis) 

All  the  output  illustrated  in  paragraph  5.  3.  3  is  also  available  for  an  error 
analysis  (FEIGN)  run.  Example  15  illustrates  the  type  of  output  obtainable 
by  input  control.  Any  combination  of  the  matrices  may  be  selected  as  output 
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by  option.  Only  the  square  roots  of  the  diagonals  of  the  covariance  matrices 
may  be  specified  as  output,  if  desired.  Printout  occurs  at  the  times  specified 
by  PRTIM. 

5.3.7  Other  Output 

There  are  three  additional  means  of  acquiring  output,  but  since  each  requires 
more  familiarity  with  the  programming  of  TRACE,  they  will  only  be  mentioned 
here. 

First,  the  FORTRAN  DEBUG  system  may  be  used.  All  binary  routines  com¬ 
piled  from  FORTRAN  source  decks  are  preceded  by  their  symbol  table,  and 
almost  all  information  of  interest  can  be  dumped  from  COMMON. 

Second,  the  source  programs  themselves  may  be  modified  and  recompiled. 

Third,  core  dumps  may  be  obtained  in  case  of  trouble  or  at  the  end  of  a  run. 

To  this  end,  each  link  contains,  as  its  first  program,  a  copy  of  SDUMP.  A 
suitable  manual  transfer  (to  a  location  dependent  on  the  version  of  FORTRAN 
in  present  use)  will  automatically  produce  an  octal  dump. 
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21606900-  6  25460.  ORAG  .015151515  FINP  INPUT  CARO  22 
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EXAMPLE  2  (continued) 
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APPENDIX 

Standard  Values  of  Constants  and  Parameters 

Included  in  this  appendix  are  lists  of  the  contents  of  the  arrays  OBJZ, 
OBJT,  OBLT,  INTEG,  C,  NUMB,  and  IFLAG.  The  first  three  contain  the 
earth  gravitational  moc  el  used;  INTEG  contains  the  numerical  integration 
parameters;  C  contains  various  constants.  The  non- zero  values  given  for 
these  arrays  comprise  a  'standard'  set  of  values  to  be  input.  The  two  arrays 
NUMB  and  IF LAG  contain  many  items  that  are  equivalent  to  input  items 
described  in  Section  5.  There  are  also  additional  items  that  the  user  may 
want  to  alter. 
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X  IS  THE  LOGICAL  TAPE  UNIT  ON  WHICH  THE 
TRAJECTORY  IS  WRITTEN. 

THE  CURRENT  TRAJECTORY  IS  DIFFERENCED  WITH  THE 
ONE  ON  TAPE  UNIT  X  AND  THE  OIFFERENCES  ARE 
PLACED  ON  TAPE  UNIT  X+l. 
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